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Abstract 

This  research  effort  develops  a  program  using  MATLAB  to  solve  the  equations  of 
motion  for  atmospheric  reentry  and  analyzes  the  validity  of  the  program  for  use  as  a  tool 
to  expeditiously  predict  reentry  profiles.  The  reentry  vehicle  is  modeled  as  a  point  mass 
with  constant  aerodynamic  properties  as  defined  by  the  user.  The  equations  of  motion  for 
reentry  are  based  on  the  two-body  problem.  The  atmosphere  is  modeled  as  a  single  layer 
exponentially  decreasing  in  density.  The  MATLAB  program  has  the  ability  to  derive  the 
initial  trajectory  conditions  from  the  position  and  velocity  relative  to  the  rotating  surface 
of  the  Earth,  the  Earth-centered  inertial  position  and  velocity,  or  the  classic  orbital 
elements.  The  program  is  compared  to  previously  established  programs  in  order  to 
validate  its  accuracy  and  numerical  stability  when  predicting  various  reentry  profiles  to 
include  sub-orbital,  super-circular  and  hyperbolic  trajectories  as  well  as  wide  ranges  of 
aerodynamic  properties. 
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DEVELOPMENT  AND  VALIDATION  OF  REENTRY  SIMULATION  USING 

MATLAB 


I.  Introduction 


Background 

Over  the  years  several  programs  have  been  developed  to  predict  the  dynamics  of 
vehicles  during  atmospheric  reentry.  The  industry  and  military  standard  is  the  Program 
to  Optimize  Simulated  Trajectories  (POST).  POST  is  a  generalized  point  mass,  discrete 
parameter  targeting  and  optimization  program.  The  program  was  originally  developed 
for  NASA  in  1970  as  a  Space  Shuttle  trajectory  optimization  program.  Since  that  time  it 
has  been  repeatedly  updated  to  improve  its  capabilities  in  the  areas  of  vehicle  modeling, 
trajectory  simulation  and  targeting  optimization.  It  features  modular  program 
construction  and  generalized  routines  so  that  various  vehicle  and  planetary  parameters 
can  be  simulated.  It  is  based  on  FORTRAN  77  and  C  programming  languages  (Brauer 
and  others,  1977:1-3). 

POST  provides  the  capability  to  optimize  point  mass  trajectories  for  both  powered 
and  unpowered  vehicles  throughout  an  entire  regime  of  ascent  to  orbit  and  maneuvers  to 
reentering  the  atmosphere  of  a  rotating,  oblate  or  spherical  planet.  The  program  integrates 
the  equations  of  motion  for  a  user-defined  model  while  solving  a  constrained 
optimization  problem  to  maximize  or  minimize  a  user  defined  function.  The  user 
supplies  a  guess  for  the  optimal  control  vector  and  POST  iteratively  solves  the  problem 
for  the  optimal  control  vector  that  maximizes  a  scalar  objective  function  subject  to  the 
user  defined  constraints  (Brauer  and  others,  1977:7,8). 
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Typical  applications  of  POST  include  a  variety  of  simulations  such  as  ascent  to 
orbit  with  optimization  of  payload,  velocity  or  initial  weight;  ascent  abort  with 
maximization  of  abort  interval;  intercontinental  ballistic  missiles  with  maximization  of 
payload;  reentry  with  optimization  of  heat  rate  or  crossrange  distance;  orbital  maneuvers 
to  minimize  fuel  consumption;  and  aircraft  performance  to  maximize  velocity,  cruise 
time  or  payload  (Brauer  and  others,  1977:13).  Among  other  projects,  POST  has  recently 
been  used  in  the  planning  for  the  Mars  Airplane  (Murray,  2001:3),  the  aerocapture 
simulation  for  the  Titan  Explorer  Mission  to  the  Saturnian  system  (Way,  2003:1)  and 
NASA’s  Reusable  Launch  Vehicle. 

A  program  called  TARGET  was  developed  as  a  simpler  and  faster  alternative  to 
POST.  TARGET  was  created  for  use  as  a  tool  in  the  systems  concept  stage  where  speed 
and  reasonable  accuracy  are  sufficient  for  planning.  Although  POST  can  model  much 
more  complicated  motion  than  TARGET,  trajectories  simulated  by  TARGET  have  been 
validated  as  virtually  identical  to  the  results  POST  generates  with  the  same  input 
parameters.  TARGET  also  generates  a  trajectory  in  1/1 0th  the  time  required  by  POST. 
As  with  POST,  TARGET  is  programmed  in  FORTRAN.  TARGET  does  show  some 
numerical  instability  when  unrealistic  final  velocities  are  targeted  and  when  extreme 
dynamics  such  as  hyperbolic  velocity  combined  with  negative  lift  are  modeled 
(Hicks, 1993:4-1). 

TARGET  solves  the  equations  of  motion  for  reentry  as  a  two-point  boundary- 
value  problem.  It  uses  a  “double-loop”  shooting  algorithm  to  optimize  the  simulated 
trajectory.  The  “double-loop”  method  attempts  to  match  the  unspecified  “free” 
conditions  to  the  boundary  conditions  identified  by  the  user  by  starting  with  an  initial 
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guess  and  iteratively  solving  the  two-point  boundary-value  problem  until  an  optimal 
solution  is  found  and  the  “free”  conditions  are  identified.  Due  to  the  fact  that  it  is  a 
targeting  program,  the  user  must  identify  the  final  radius,  longitude,  latitude,  and  velocity 
for  the  reentry  vehicle.  Additionally,  the  user  must  supply  the  initial  radius,  flight-path 
angle  relative  to  the  plane  perpendicular  to  the  radius,  and  the  orbital  inclination.  The 
program  is  able  to  generate  a  reentry  profile  to  include  variations  in  velocity,  altitude, 
flight-path  angle,  latitude,  longitude,  heading,  Mach  number,  and  rate  of  heating  (Hicks, 
1993 :  A- 1). 

A  third  program  of  similar  nature  to  the  objective  of  this  research  is  called 
IMPULSE.  The  National  Air  and  Space  Intelligence  Agency  (NASIC)  developed 
IMPULSE  as  a  MATLAB  based  means  to  simulate  ballistic  missile  flight  paths. 
IMPULSE  can  model  the  entire  trajectory  of  missiles  from  boost  through  impact.  The 
unclassified  version  of  the  program  contains  examples  of  missiles  and  RVs  that  are  not 
classified  and  is  preloaded  with  the  geographical  data  of  multiple  launch  and  impact  sites 
across  the  globe.  Multiple  reentry  vehicles  (RVs)  can  be  simulated.  It  also  includes  the 
ability  to  model  events  during  the  missile  flight  that  change  the  mass  or  aerodynamics  of 
the  RV,  such  as  booster  separation,  shroud  ejection,  or  change  in  angle  of  attack 
(IMPULSE  2005). 

The  focus  of  this  research  is  the  development  and  validation  of  a  reentry  simulator 
using  the  same  language  (MATLAB)  as  IMPULSE.  Whereas  IMPULSE  is  a  complex 
application  that  considers  the  entire  trajectory  of  a  ballistic  missile  with  only  a  rough 
analysis  of  the  reentry  portion,  this  research  will  focus  on  a  quick  solution  for  the 
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simulation  of  a  wide  variety  of  entry  trajectories  with  much  smaller  integration  steps  than 
IMPULSE  employs. 

MATLAB  (MATrix  LABoratory)  is  a  high-level  technical  computing  language 
for  algorithm  development,  data  visualization,  and  data  analysis  intended  for  use  with 
numerical  computations.  The  software  has  a  wide  range  of  applications,  including  signal 
and  image  processing,  communications,  control  design,  test  and  measurement,  financial 
modeling  and  analysis,  and  computational  biology.  The  original  goal  was  to  enable 
scientists  and  engineers  to  use  matrix-based  techniques  to  solve  problems,  without  having 
to  write  programs  in  traditional  languages  like  C  and  FORTRAN.  Over  the  years 
additional  capabilities  have  been  added,  in  particular  outstanding  graphics  support. 
MATLAB  is  rapidly  displacing  other  languages  as  the  software  of  choice  for  algorithm 
programming  in  education  and  industry  and  is  now  the  preferred  tool  for  many  agencies 
of  the  U.  S.  Air  Force  and  the  Department  of  Defense  (General  2006). 

This  paper  describes  the  development  and  validation  of  the  trajectory  simulation 
program.  Chapter  II  begins  with  a  description  of  the  reference  frames  used  to  describe 
the  position  and  direction  of  travel  of  the  reentry  vehicle,  followed  by  development  of  the 
equations  of  motion.  Chapter  III  shows  how  the  initial  conditions  were  established  for 
the  equations  of  motion.  Chapter  IV  presents  a  comparison  of  the  MATLAB  program 
developed  in  this  research  with  established  software  programs  and  analytical  solutions. 
A  description  of  the  program  operation  is  also  included.  The  final  chapter  (V)  states  the 
conclusions  drawn  from  the  research  and  provides  recommendations  for  improvements  to 
the  software  that  may  be  undertaken  in  the  future.  Finally,  the  Appendices  include  a 
flowchart  of  the  alogorithm,  descriptions  of  the  MATLAB  code  and  the  code  itself. 
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II.  Equations  of  Motion  for  Reentry 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  explain  the  development  of  the  equations  of 
motion  for  atmospheric  reentry.  First,  the  reference  frames  used  in  this  research  will  be 
established.  Next,  the  equations  of  motion  will  be  derived  for  a  point  mass  flying  into  the 
atmosphere  of  a  rotating  planet.  All  of  the  equations  and  figures  in  this  chapter  are 
summarized  from  Hicks  (2006:33-58). 

Reference  Frames 

The  geocentric-equatorial  system  is  an  inertial  reference  frame  with  its  origin  at 
the  center  of  the  planet.  In  this  system  the  x-y  plane  is  the  equatorial  plane  with  the  x- 
axis  pointing  at  the  vernal  equinox.  The  z-axis  passes  through  the  north  pole  of  the 
planet.  The  y-axis  is  positioned  such  that  e_  =etxev.  In  the  case  of  the  Earth  this  is 

commonly  referred  to  as  the  Earth-Centered  Inertial  (ECI)  reference  frame  and  will  be 
referred  to  as  such  in  this  research. 

In  the  planet-fixed  coordinate  frame  the  system  spins  about  the  z-axis  with  the 
rotation  of  the  planet.  In  the  case  of  the  earth  this  reference  frame  will  be  referred  to  as 
the  Earth-Centered  Fixed  (ECF)  coordinate  system.  In  the  ECF  system  the  x-axis  points 
through  the  Greenwich  Meridian  (0°  longitude).  The  relation  between  the  ECI  and  ECF 
coordinate  systems  is  illustrated  in  Figure  (1)  where  the  ECF  system  is  shown  as 

OXiYjZi  the  planet  rotates  with  a  constant  angular  velocity  a> and  At  is  the  sidereal 
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time.  Sidereal  time  is  defined  as  the  time  elapsed  since  the  two  systems  were  last 
aligned. 


Cz,  Cz, 


Figure  1.  ECF  and  ECI  Coordinate  Frames 


A  coordinate  transformation  from  the  ECI  frame  to  the  ECF  frame  is  given  by: 


cos(<z>Af) 

-sin(ft)At) 

0 


sin(<z>Af)  0 
cos(<z>At)  0 
0  1 


[*] 


(i) 


A  third  system  also  originates  at  the  planet’s  center  but  defines  the  x-axis  as 
pointing  along  the  position  vector  (r  )  of  the  vehicle  of  interest.  As  shown  in  Figure  (2), 
this  frame  (OX2Y2X2)  can  be  created  by  a  rotation  9  about  the  ECF  z-axis  followed  by  a 
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rotation  -(f)  about  the  y-axis.  Therefore,  0  represents  the  longitude  and  (f)  represents 


the  geocentric  latitude. 


Figure  2.  ECF  and  Vehicle-Pointing  Coordinate  Frames 


Performing  the  two  rotations  yields  the  following  relationship  between  the  ECF 
frame  (OX1Y1Z1)  and  the  vehicle-pointing  frame  (OX2Y2Z2): 


cos  (f)  cos  6 
-sin  9 
-sin^cos# 


cos  (j)  sin  0 
cos  0 

-  sin  (f>  sin  9 


sin  (j) 

0 


COS  (j) 


(2) 


Figure  (3)  introduces  two  more  angles  necessary  to  help  describe  the  trajectory  for 
reentry.  The  flight-path  angle,  Y,  is  defined  as  the  angle  between  the  velocity  vector,  V  , 
and  the  plane  orthogonal  to  r  (the  local  horizontal  plane).  Y  is  positive  when  V  is  above 
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the  local  horizontal  plane.  The  heading  angle,  t// ,  is  defined  as  the  angle  between  the 


local  parallel  of  latitude  and  the  projection  of  V  on  the  local  horizontal  plane,  i//  is 
positive  in  the  right-handed  direction  about  the  X2-axis. 


Figure  3.  Flight  path  and  heading  angles  in  relation  to  ECF  and  Vehicle-Pointing  Systems 

A  useful  relationship  to  explain  at  this  point  is  one  that  relates  angular  motion 
(Q)  between  the  ECF  system  and  the  vehicle -pointing  system.  The  vehicle-pointing 
system  was  created  by  a  rotation  0  about  the  zi-axis  followed  by  a  -(f)  rotation  about  the 
y2-axis.  Therefore,  the  angular  rate  can  be  written  as: 
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(3) 


Equation  (2)  can  then  be  used  to  write  this  expression  entirely  in  terms  of  the  vehicle¬ 
pointing  system: 


Q  = 


A 


dd  .  A 
—  sin^ 

v  dt  j 


d(j) 

dt 


e_v2  + 


A 


de 

- COS0 

v  dt  j 


(4) 


Equations  of  Motion 

By  examining  Figure  (3)  one  can  derive  the  equation  for  the  velocity,  denoted 
here  as  RV  .  The  “R”  superscript  denotes  that  the  velocity  is  relative  to  the  rotating  ECF 
frame: 

RV  =  (^Fsiny)  +  (fiFcos;rcos^j  eVi  +  ( RV  cosy  sin  i//J  e^  (5) 

A  derivative  in  ECF  frame  is  related  to  a  derivative  in  vehicle-pointing  reference 
frame  by  the  following  expression 


'dr 

dt 


2  dr 
dt 


+  co 


2/1 


x  r 


(6) 


i 


where  —  is  the  velocity  relative  to  the  ECF  system  and  com  is  the  angular  rate 
dt 

established  in  Equation  (4)  as  Q  .  Therefore,  Equation  (5)  can  be  rewritten  as: 
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2 


(V) 


- 1-  Q  x  r 

dt 


After  completing  the  cross-product  and  writing  in  terms  of  the  vehicle-pointing  system. 
Equation  (7)  becomes: 

RV  =  r  eX2  +(r0 cos  <p)  +  (r^)  eZj  (8) 

Equations  (5)  and  (8)  are  both  solutions  for  the  relative  velocity  written  in  tenns  of  the 
vehicle-pointing  system.  Therefore,  each  component  of  Equation  (5)  must  be  equivalent 
with  its  corresponding  component  in  Equation  (8): 


r=RV  sin  y 

(9) 

x  RV  cosy  cos  y/ 
u  = 

rcostf) 

(10) 

^  RV  cosy  sin  y/ 

(ID 

r 


These  are  known  as  the  three  kinematic  equations  of  motion,  which  define  the  vehicle’s 
position  relative  to  the  rotating  planet. 

These  equations  only  give  half  of  the  information  needed  to  fully  describe  the 
motion  of  the  vehicle  during  reentry.  Three  more  expressions  that  account  for  the 
changes  in  velocity  and  direction  due  to  the  forces  of  gravity  and  aerodynamics  are 
needed. 
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For  the  purposes  of  this  research  the  reentry  vehicle  is  treated  as  a  thrustless  point 
mass.  The  total  force  acting  on  the  point  mass  can  then  be  defined  as  the  summation  of 
the  force  from  gravity  and  the  aerodynamic  forces: 

F  =  A  +  mg  (12) 

Another  definition  of  the  total  force  comes  from  Newton’s  Second  Law  where  it  is 
assumed  the  mass  is  constant  and  the  reference  frame  is  inertial 

F  =m  -  (13) 

dt 

where  the  “I”  in  the  superscript  shows  that  this  is  a  derivative  of  the  velocity  in  the 
inertial  frame.  In  order  to  express  the  derivative  of  the  velocity  in  terms  of  the  vehicle¬ 
pointing  system  it  is  necessary  to  apply  the  same  technique  used  to  establish  the 
relationship  of  a  derivative  in  a  fixed  frame  to  a  derivative  in  a  rotating  frame  used  for 
Equation  (6) 


/  -  R  - 

dr  dr  —  - 

—  = - hcoxr 

dt  dt 


(14) 


and  taking  a  second  derivative  to  find  the  derivative  of  the  inertial  velocity 


‘  d2r  _  !dV 
dt2  dt 


*  d 

(  R 

X 

(  R 

dr  —  - 

— 

dr  —  - 

— 

- b  cox  r 

+  &>x 

- b  cox  r 

dt 

dt 

dt 

V 

J 

V 

J 

(15) 
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where  co  is  the  angular  rate  between  the  ECI  and  ECF  frames.  For  the  purpose  of  this 


— -  Rdr 

research,  this  angular  rate  can  be  treated  as  a  constant  called  cor]j .  -  is  the  velocity 

dt 

Rd2r 

relative  to  the  rotating  surface  of  the  planet  and  — —  is  the  acceleration  relative  to  the 

dt~ 

rotating  surface  of  the  planet.  Equation  (15)  simplifies  to: 


dV  dRV 


dt  dt 


+  2 


(®ex^)  +  ®®x(fflexr) 


(16) 


dV  F 

Replacing  - with  —  from  Equation  (13)  results  in: 

dt  in 


m 


Rd_  V 
dt 


( - -  R-~\ 

“ - ►  / - ► 

oe  x  V  -  m 

®®x(0ex/-  J 

(17) 


As  defined  in  the  ECF  system,  co®  is  about  the  z-axis  only.  It  can  be  converted  to  the 
vehicle-pointing  system  by  applying  Equation  (2): 

rue  =(oesin^)ex2 +(®ecos^)e_-2  (18) 

- *  ^  ^ 

Equations  (5)  and  (18)  can  then  be  used  to  find  the  cross-product  com  x  V  in  Equation 
(17): 


- "  ^  - -  /  R  \  ^ 

r%x  V  =  -y  Vcom  cos  (j)  cosy  cos  ^)e.t2 

+  RVa>m  (cos  </> sin  y- sin  ^cosysin^e^  (19) 

+  ( RVcve  sin  </>  cos  y  cos  y/  )  eZ2 
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Similarly,  recognizing  that  r  =  rc:<2  the  cross  product  com  x  [  x  rj  can  be  solved: 

x(®e  xr)  =  (ro^  cos2  ^)e.v2  +(ra>®  sin^cos^)eZ2  (20) 

Recall  that  the  forces  ( F )  acting  on  the  vehicle  are  a  summation  of  the 
aerodynamic  forces  (A)  and  the  force  due  to  gravity  ( mg ).  The  aerodynamic  forces  can 

be  further  broken  down  into  a  drag  force  ( D  )  and  a  lift  force  ( L ).  Drag  is  defined  as 
acting  in  the  direction  opposite  velocity,  and  can  therefore  be  described  from  Figure  (3) 
in  the  vehicle-pointing  system  by  geometric  relation  as: 

D  =  -(Dsiny)e.v2  -(Dcosy  cos^)ey2  -  (D  cosy  sin  ^)e_-2  (21) 

In  order  to  define  the  lift  force  a  new  variable  describing  the  bank  angle,  cr ,  must 

-  R  — 

be  introduced,  cr  is  defined  as  the  orientation  of  the  lift  vector  relative  to  the  ( r,  V  ) 
plane  as  shown  in  Figure  (4). 
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Recognizing  that  the  lift  vector  is  also  perpendicular  to  the  drag  vector,  lift  can 
then  be  described  in  the  vehicle-pointing  system  as  a  function  of  bank  angle,  flight-path 
angle  and  heading: 

L  =  L  (cos  <j  cos  y )  e*2  +  L(-  cos  a  sin  y  cos  yr  -  sin  a  sin  y/)eyi 
+L(-  cos  cr  sin  y  sin  y/  +  sin  cr  cos  ^)e_-2 


The  final  component  of  the  force  vector  ( F  )  is  the  gravity  term  mg  .  Since  the 
bodies  in  the  system  are  being  treated  as  point  masses  for  the  purpose  of  this  research, 
gravity  always  acts  along  the  radius  vector: 


mg  =  (-fflg)e. 


(23) 


Now  all  of  the  terms  on  the  right  hand  side  of  Equation  (17)  have  been  described 
in  the  vehicle-pointing  system.  The  equation  can  be  rewritten  as: 


m 


Rd[RV ) 
dt 


+  L( cos  cr  cos  y)  -  (D  sin  y )  -  mg 


-2m 


1  e„ 


-(rV(d9  cos  (j>  cosy  cos  y/^j 
- m\^-{r(0^cos2 

( +L(- cos  cr  sin  y  cos  y/ -s\n  a  smy/)  + D(- cos  y  cos  i//)j 
|-2m  rVco9  ( cos  (f>  sin  y  -  sin  (j>  cosy  sin  y/) 

\ +L(- cos  cr  sin  y  sin  y/  +  sin  cr  cos  y/)  +  D  (-  cos  y  sin  y/ )  | 

1 -2w[(  rVg>9  sin  (j>  cosy  cos  ^)J-w[rcy^  sin^cos^  ] 


(24) 
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It  is  now  necessary  to  equate  the  derivative  of  the  velocity  relative  to  the  ECF 


frame 


Rd[Rv) 

dt 

V  J 


with  the  derivative  of  the  velocity  in  the  vehicle-pointing  system.  The 


same  technique  used  in  Equations  (6)  and  (14)  is  applied  here: 


v(v)  2j(y) 


dt 


dt 


+  Qx*F 


(25) 


where  Q  is  the  angular  rate  defined  in  Equation  (3). 
Equation  (5)  to  yield: 


dt 


is  easily  obtained  from 


2d(Rf ) 

— ^ — -  =  (RVsiny  +  RVycosy"j  eVj 

+  cosf  cos^-  RVysinycosy/~  RVyr  cosy  sin^)  e 
+  ( RV cosy siny/ -  RVy  siny  siny/  +  RV y/ cos y cos^)  e 


^2 


(26) 


The  cross-product  QxB  V  from  Equations  (4)  and  (5)  is: 


f  • 


QxRV=RV 


+RV 


-  (j)  cos  y  sin  y/  -  9  cos  (j)  cos  y  cos  y/ 


0  cos  (j)  sin  y  -  6  sin  (f>  cos  y  sin  yr 


v 

+  Rv[  6  sin  (j)  cosy  cosy/  +  (j)siny 


•'yi 


(27) 


Replacing  6  and  (j)  with  the  expressions  defined  in  the  kinematic  Equations  (10)  and 
(11),  and  combining  Equations  (25)  -  (27)  gives: 
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Rd{Rv) 

dt 


.  V2C0SV'' 


Vsin/+  V  y  cosy- - 


?E  cos  7  cos  y/  -  RVy  sin  y  cos  y/-RVy/  cos  7 sin  y/ 

RV2 

H - cos  y  cos  ^  ( sin  y  -  tan  ^  cos  y  sin  ^ ) 

r 

v  cosy  siny/  -  RV  y  smy  smy/  +  RVy/ cosy  cosy/ 


>’2 


Ry2 


cosy  (cosy  cos2  ^/■tan^  +  sinysin^/J 


(28) 


Equations  (24)  and  (28)  are  both  solutions  for  the  derivative  of  the  velocity  in  the 
ECF  frame  expressed  in  vehicle -pointing  coordinates,  and  can  therefore  be  equated 
component-by-component  to  generate  three  coupled  differential  equations.  After  a  bit  of 

R  .  . 

work,  they  can  be  solved  for  V,  y  and  y/ : 


RV  = 


D 

m 


g  sin  y  +  rcoh  cos  (f>  ( cos  (ft  sin  y  -  sin  (ft  sin  y/  cos  y ) 


(29) 


RVy 


L  RV2 

=  H — cos  <J  —  g  cos  y  h - cosy  +  2  cos  (ft  cos  y/ 

m  r 

+rcol)  cos^(cos0cosy  +  sin0sin^siny) 


(30) 


R  .  L  sin  (7 


Vy/  = 


Ry2 


m  cos  y 


cos  y  cos  y/  tan  (ft +  2  RV co e  (sin  y/  cos  (ft  tan  y  -  sin  (ft) 


(31) 


roh 


cos  y 


-  sin  ^  cos  ^  cos  y/ 
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These  are  the  three  differential  force  equations.  When  integrated  simultaneously 
with  the  three  kinematic  equations,  the  location  and  movement  of  a  point  mass  in  a  two- 
body  system  under  the  influence  of  aerodynamic  forces  can  be  completely  described  by 
the  six  independent  parameters  ( r ,  9,  <f> ,  RV  ,  y  and  i//  ).  Note  that,  just  as  the  velocity 
(  RV  )  is  relative  to  the  rotating  atmosphere,  so  is  flight-path  angle  ( y )  and  heading  ( i//  ). 

The  focus  of  this  research  is  to  set  up  a  MATLAB  program  to  integrate  these 
equations  and  to  test  their  validity  and  stability  at  various  trajectories. 
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III.  Problem  Setup 


Overview 

This  chapter  provides  details  of  the  practical  setup  of  the  problem  for  a  point- 
mass  entering  the  Earth’s  atmosphere.  It  includes  definitions  of  the  terms  and  the 
methods  used  to  derive  the  initial  values  for  the  equations  of  motions  from  ECI 
coordinate  and  also  from  classic  orbital  elements. 

Defining  Terms 

As  stated  earlier  in  the  derivation  of  Equation  (16),  the  angular  rate  of  the  Earth’s 
rotation  is  assumed  to  be  constant.  The  value  for  com  can  then  be  defined  (WGS  84 
1998:14): 


0)@  =7.2921 15855228083A- 5  rad/s 


(32) 


The  gravity  term  (g),  on  the  other  hand,  is  not  a  constant  but  varies  as  a  function 
of  the  radius  (Regan,  1984:23) 


g  =  g  0 


(r  V 

'o 

V  r) 


(33) 


where  the  force  of  gravity  and  the  radius  of  the  Earth  at  the  equator  are  the  constants  g0 
and  r0,  respectively.  The  WGS-84  values  for  g0  (Department  2000:3-7)  andr0 
(Department  2000:3-2)  are: 

g0  =9.7803253359  m/s2  (34) 
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(35) 


r0  =  6378137  m 

The  values  for  lift  (L)  and  drag  (D)  are  approximated  by  the  following 
expressions  (Hicks,  1993:3-5) 


D  = 


pSCD  Rv2 
2 


(36) 


L  = 


pSCL  Rv_ 2 
2 


(37) 


where  S  is  the  reference  surface  area  of  the  vehicle,  CD  and  C,  are  the  coefficients  of 
lift  and  drag,  and  p  is  the  density  of  the  atmosphere.  For  this  research,  the  assumption  is 
made  that  the  reference  surface  area  and  aerodynamic  coefficients  are  constant  and 
known  for  the  reentry  vehicle.  The  density  of  the  atmosphere  varies  with  altitude.  For 
simplification,  this  study  uses  a  single  layer  exponential  atmosphere  described  by 
(Duncan,  1962:124) 


p  =  P^ph 

(38) 

p0  =  1.54793  kg/m3 

(39) 

fi  = - - - 

(40) 

6935  m 


where  p0  is  the  density  of  the  atmosphere  at  sea-level,  ft  is  the  inverse  of  the  scale 

height  for  the  atmosphere  and  h  is  the  altitude  in  question.  The  surface  density  and  scale 
height  are  variables  by  temperature  and  other  factors,  but  average  values  were  chosen  for 
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simplicity  (Hicks,  1993:4-9).  This  exponential  model  shows  excellent  agreement  with 
the  1976  Standard  Atmosphere  to  an  altitude  of  approximately  120  km,  beyond  which 
these  values  are  far  more  difficult  to  model  due  to  large  fluctuations  (Regan  1993:39). 

Determining  the  altitude  is  easy  if  the  Earth  is  treated  as  a  sphere.  It  would 
simply  be  the  difference  between  the  radius  (r)  to  the  vehicle  and  the  constant  radius  of 

the  Earth.  But  the  surface  of  the  Earth  is  actually  an  uneven  geoid,  which  is  better 
modeled  by  an  ellipsoid  whose  ellipticity  ( e )  is  defined  as  (Department,  2000:3-7) 

e  =  l~—  (41) 

>o 

rp  =6356752.3 142  m  (42) 

where  rp  is  the  polar  radius  of  the  Earth.  The  altitude  of  the  vehicle  ( h ),  shown  in 
Figure  (5),  can  then  be  approximated  if  the  assumption  is  made  that 

r-rE=h  =  l  (43) 

where  rE  is  the  radius  of  the  Earth  as  a  function  of  geocentric  latitude.  This  assumption 

is  valid  for  the  Earth  because  its  ellipticity  is  very  small  (« 1/298)  and  yields  a  result 
accurate  to  within  2  m  for  an  altitude  of  300  km  (Regan,  1984:27).  Now,  a  very  good 
approximation  for  rE  is 


cos  (f) 


(44) 
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where  xa  is  the  distance  shown  in  Figure  (5)  and  described  by  (Geocentric  2005): 


Figure  5.  Reference  Ellipsoid 
Establishing  the  Initial  Conditions 

For  this  study,  the  differential  equations  of  motion  were  solved  as  initial  value 
problems  using  MATLAB.  In  order  to  set  the  initial  conditions,  it  is  essential  to  define 
the  values  needed  for  the  differential  equations.  In  the  unlikely  case  the  user  knows  these 
particular  parameters,  they  can  be  entered  directly  as  the  initial  conditions.  If  other  more 
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common  values  are  known  for  the  reentry  vehicle,  conversion  to  the  initial  values 
required  for  the  integration  is  necessary. 

As  a  means  to  provide  flexibility  for  the  user,  the  program  was  written  to  allow 
for  inputs  from  three  different  coordinate  systems.  First,  the  user  may  enter  the  initial 
values  for  the  differential  equations  directly  in  terms  of  the  vehicle-pointing  system 
coordinates  relative  to  the  ECF  frame.  Second,  the  values  for  position  and  velocity  may 
be  defined  in  terms  of  the  Cartesian  ECI  frame.  Finally,  the  six  classic  orbital  elements 
can  be  used  to  define  the  initial  reentry  conditions.  In  all  cases  the  user  must  also  define 
the  constants  for  the  reentry  vehicle’s  mass,  coefficients  of  lift  and  drag,  the  effective 
aerodynamic  surface  area  and  the  constant  vehicle  bank  angle. 

Finding  Initial  Values  from  Inertial  Position  and  Velocity 

If  the  initial  conditions  are  given  in  the  ECI  frame,  a  method  to  convert  them  to 
the  values  desired  follows.  The  first  step  is  to  calculate  the  orbital  momentum  vector 

(H)  (Weisel,  1997:62): 


H  =  r  x  v  (46) 

The  inertial  flight-path  angle  relative  to  the  local  horizontal  plane  can  then  be 
found  (Hicks,  2006:17) 


cos  Iy  =  —  (47) 

rv 

which  is  ambiguous  in  terms  of  the  sign  of  the  inertial  flight-path  angle.  The  flight-path 
angle  is  defined  as  negative  when  it  is  Earthward  of  the  local  horizontal  plane;  therefore, 
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if  the  vehicle  is  moving  away  from  the  origin  of  the  ECI  frame  the  flight-path  angle  is 


positive.  In  other  words,  if  r  •  v  >  0 ,  the  flight-path  angle  is  positive;  otherwise,  it  is 
negative. 

The  initial  latitude  ((/))  is  simple  to  find  due  to  the  fact  that  the  xy-plane  in  the 
ECI  frame  corresponds  to  the  Earth’s  equator 


sin^  =  —  (48) 

r 

where  xyz  are  the  components  of  r.  Since  -90°  <(f)<  90° ,  the  above  equation  detennines 
(f>  without  ambiguity. 

The  longitude  (6)  is  the  angle  between  the  position  vector  projected  into  the  xy- 
plane  and  the  positive  x  direction: 


sin  6  =  —j=^=  (49) 

V*  +y 

This  will  only  yield  results  of  ±90° ,  so  a  correction  can  be  made  by  subtracting  the  result 
for  0  from  n  when  the  x  component  of  the  initial  position  is  negative.  Equation  (49)  is 
only  valid  when  the  ECI  and  ECF  coordinate  systems  are  aligned.  Otherwise,  a 
correction  must  be  made  for  the  angular  position  of  the  ECF  frame  relative  to  the  ECI 
frame: 


6  =  sin 


A 


y 


{Jx2+y2  ) 


—  Ato) 


(50) 
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The  inclination  of  the  trajectory  is  needed  in  order  to  calculate  the  heading  of  the 


reentry  vehicle.  Inclination  is  the  angle  between  the  angular  momentum  vector  ( H  )  and 
the  z-axis,  given  by  (Weisel,  1997:62): 


COS  l  = - 

H 


(51) 


The  inertial  heading  (V)  is  the  angle  between  the  projection  of  the  inertial 
velocity  onto  the  local  horizontal  plane  and  the  latitude.  The  inertial  heading  can  be 
related  to  the  inclination  and  latitude  (Longusky  and  Vinh,  1980:8): 


/  cos  i 

cos  V  = -  (52) 

COStf) 

This  result  has  quadrant  ambiguity  because  the  heading  can  vary  a  full  360°.  In  order  to 
resolve  this  problem,  it  is  necessary  to  find  the  right  ascension  of  the  ascending  node. 
First  the  line  of  nodes  n  can  be  found  as  the  vector  perpindicular  to  the  angular- 
momentum  vector  H  and  the  z-axis  (Weisel,  1997:62): 


n  = 


Hxe: 

Hxe: 


(53) 


The  line  of  nodes  lies  in  the  equatorial  plane  and  has  the  form: 


n  =  cx  cos  Q  +  e  v  sin  Q 


(54) 
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The  right  ascension  of  the  ascending  node  ( Q  )  is  now  known  without  ambiguity.  The 
inertial  heading  is  always  positive  when  the  object  in  question  is  within  ±90"  of  the 
ascending  node.  Therefore,  when  the  difference  between  0 ,  found  in  Equation  (49),  and 
Q  is  less  than  |±90°| ,  the  result  is  0°  <  V  <  180" ,  otherwise  180°  <  'y/  <  360° . 

The  next  task  is  to  find  the  heading  relative  to  the  rotating  atmosphere.  For 
simplification,  it  is  a  good  idea  to  define  the  local  horizontal  component  of  the  inertial 
velocity  as  Vh 


Vh  =  V  cos  ‘y  (55) 

and  the  velocity  of  the  rotating  atmosphere  ( Ve)  at  the  position  of  the  vehicle  relative  to 
the  center  of  the  Earth  (Longusky  and-  Vinh,  1980:8): 

Ve=rco9  cost/)  (56) 


The  horizontal  component  of  the  relative  velocity  is  then  Vh-Ve.  The  relative 
heading  y/  is  the  angle  between  the  horizontal  component  of  the  relative  velocity  and  the 
local  parallel  of  latitude  as  shown  in  Figure  (6).  y/  can  be  found  from  the  law  of  cosines 


cos^-  = 


K  cos  >-Fg 

^v,1  +  vt1-2v,rt<m‘w 


(57) 


where  the  quadrant  ambiguity  of  y/  is  solved  in  the  same  manner  as  for  V  . 
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Plane 


Local  parallel  of 
latitude 


Figure  6.  Relationship  Between  Heading  and  Velocity 

R  — 

The  next  value  that  needs  to  be  defined  is  the  velocity  (  V )  relative  to  the 
rotating  atmosphere  shown  in  Figure  (7).  The  relative  velocity  is  described  as  V  —  Ve  and 
once  again  the  law  of  cosines  in  invoked  to  come  up  with  a  relationship  for  the  magnitude 
of  RV : 


V  =  Vl2+L2-2LL,cos  V  (58) 

As  described  in  Chapter  1  and  illustrated  in  Figure  (7),  the  relative  flight  path 
angle,  y ,  is  defined  as  the  angle  between  the  relative  velocity  and  the  local  horizontal 
plane.  It  is  detennined  by  the  relationship: 
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(60) 


sin  /  = 


V  sin  ‘y 
RV 


Since  -90°  <  /  <  90° ,  Equation  (60)  detennines  /  without  ambiguity. 


Figure  7.  Visualizing  Relative  Velocity  and  Flight-path  Angle 


All  of  the  initial  values  for  the  equations  of  motion  (r,<j),0,ii/,  RV,yj  have 


now 


been  found  and  can  be  integrated  to  solve  for  the  entry  trajectory. 


Finding  Initial  Values  from  Classic  Orbital  Elements 

As  mentioned  earlier  in  this  chapter,  the  initial  values  for  the  equations  of  motion 
can  also  be  derived  from  the  six  classic  orbital  elements.  The  semimajor  axis  (a)  defines 
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the  size  of  the  orbit  and  its  period.  The  eccentricity  (e)  determines  the  orbit’s  shape.  The 
right  ascension  of  the  ascending  node  (Q)  was  addressed  in  the  previous  section  and 

identifies  the  points  where  the  orbital  plane  intersects  the  equator,  measured  as  the  angle 
from  the  vernal  equinox  to  the  point  where  the  orbit  crosses  the  equator  from  south  to 
north.  The  orbital  inclination  (I)  was  also  discussed  in  the  previous  section  and  is  defined 
as  angle  between  the  orbital  plane  and  the  equator  at  the  ascending  node.  The  argument 
of  perigee  (<u)  is  the  angle  from  the  ascending  node  to  the  perigee  point.  Finally,  the 

true  anomaly  (v)  is  described  as  the  position  of  vehicle  at  a  particular  instant  in  time  with 
respect  to  the  argument  of  perigee  (Weisel,  1997:60).  Most  of  these  elements  are 
illustrated  in  Figure  (8). 


Qz 
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As  the  first  step,  the  scalar  radius  can  be  directly  calculated  from  the  orbital 


elements,  as  shown  by  Griffin  and  French  (2004:106): 

a(l-e2) 

r  =  — - - - 

1  +  ecosv 


(61) 


Equation  (61)  accounts  for  all  possible  orbits  except  parabolic.  The  radius  vector  can 
then  be  written  in  the  perifocal  system,  where  the  center  of  the  Earth  is  the  origin,  ep 
points  towards  perigee  and  eQ  is  in  the  orbital  plane  perpendicular  to  ep  as  illustrated  in 
Figure  (9): 

r  =  (rcosv)ep  +(rsinv)ee  (62) 


Figure  9.  Perifocal  Coordinate  System 
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The  first  step  in  obtaining  the  velocity  vector  is  calculating  the  semi-latus  rectum 


of  the  orbit  (p)  (Griffin  and  French,  2004: 107): 

p  =  a  ( 1  -  e1 ) 


The  velocity  vector  can  then  be  found  by  (Hicks,  2006:29) 


V  = 


(-sinv)ep 


+  (e  +  cosv)e!? 


where  p  is  the  gravitational  parameter  of  the  Earth  (Weisel,  1997:323): 


(63) 


(64) 


p  =  3.98601x  1014m3  / s2  (65) 

The  values  can  now  be  found  in  the  inertial  frame  by  a  rotation.  The  substitutions 
for  the  perifocal  unit  vectors  into  the  ECI  frame  are  given  by  (Griffin  and  French, 
2004:122): 


ep  =  (cosQcos<y-sinQsin<ucosz)e.v 

+  ( sin  Q  cos  co  +  cos  Q  sin  co  cos  i )  e  v  +  ( sin  co  sin  i )  e_- 

eg  =  (- cos  Q  sin  co  -  sin  O  cos  co  cos  i)ex 

+  (-sinQsin<u  +  cosQcos(ycosz)  e  v  +  (cos  co  sin  i )  ez 


(66) 


(67) 


Now  that  the  radius  and  velocity  are  known  in  inertial  coordinates,  they  can  be 
converted  to  the  initial  values  for  the  equations  of  motion  as  described  in  the  previous 
section. 
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IV.  Program  Operation 


Overview 

This  chapter  is  a  brief  description  of  the  operation  of  the  MATLAB  reentry 
software  developed  in  this  research.  First,  a  description  of  the  user  operations  is 
presented  along  with  graphical  illustrations  of  the  interfaces.  Next,  the  software  system 
processes  are  outlined  and  presented  in  the  form  of  a  flow  chart. 

Operation 

To  begin  operation  of  the  reentry  software,  the  user  must  ensure  that  the  folder 
containing  the  program  reentry.m  is  opened  in  MATLAB’s  Current  Directory.  Typing 
“reentry”  at  the  command  line  will  then  initiate  the  program  by  bringing  up  a  graphical 
user  interface  (GUI)  titled  “Reentry  Program  Initiation”  as  shown  in  Figure  (10).  This 
GUI  allows  the  operator  to  select  the  desired  coordinate  system  of  the  initial  conditions  as 


Figure  10.  GUI  for  Reentry  Problem  Initiation 
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well  as  the  preferred  shape  of  the  atmosphere  from  drop-down  menus  listing  the  options. 
Initial  conditions  can  be  in  the  vehicle  pointing  system,  the  ECI  system,  or  described  by 
the  classic  orbital  elements.  The  atmosphere  may  be  defined  as  spherical  or  conforming 
to  the  WGS-84  ellipsoidal  approximation  for  the  Earth’s  shape. 

Clicking  “Run”  on  the  initiation  GUI  will  bring  up  a  second  GUI,  as  shown  in 
Figure  (11),  to  allow  the  user  to  input  the  vehicle  parameters  and  initial  conditions  of 
reentry.  The  operator  simply  needs  to  include  a  value  in  each  cell  corresponding  to  each 
particular  parameter.  Clicking  “Run”  initiates  the  simulation. 


Figure  11.  GUI  for  Vehicle  Parameters  and  Initial  Conditions 
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The  software  then  automatically  generates  a  plot,  like  that  shown  in  Figure  (12), 
for  each  of  the  parameters  solved  for  in  the  six  equations  of  motion;  altitude,  velocity, 
flight-path  angle,  geocentric  latitude,  longitude,  and  heading.  It  also  generates  a  plot 
showing  the  acceleration  experienced  by  the  vehicle. 


Figure  13.  Sample  Plot  of  Results 

Additionally,  the  results  are  tabulated  and  saved  into  a  matrix  named  Y.  The  Y 
matrix  can  be  selected  from  the  MATLAB  Workspace  and  displayed  in  the  Array  Editor. 
Y  has  six  columns  that  are  values  for  radius  (m),  relative  velocity  (m/s),  relative  flight- 


33 


path  angle  (rads),  longitude  (rads),  latitude  (rads),  and  relative  heading  (rads), 
respectively.  Each  row  of  Y  corresponds  to  a  time  from  the  vector  t,  which  is  also  found 
in  the  Workspace,  t  is  composed  of  the  time,  in  seconds,  for  each  step  of  the  integration, 
starting  at  zero  for  the  initial  conditions. 

The  vehicle  parameters  and  initial  conditions  are  automatically  displayed  as  a 
matrix  call  “user  inputs”  in  the  MATLAB  Array  Editor  along  with  a  matrix  of  their 
description  called  “parameters”.  The  initial  conditions  displayed  in  “user  inputs”  can  be 
modified  directly  in  the  Array  Editor  and  the  entire  simulation  can  be  run  again  by  typing 
“reentry”  at  the  MATLAB  Command  Line.  Subsequent  initiations  of  the  program  in  this 
manner  will  not  bring  up  the  GUIs  described  earlier.  In  order  to  change  the  coordinate 
system  or  atmosphere  shape,  the  MATLAB  Workspace  needs  to  be  cleared  and  the 
program  re-initiated. 

Software  System  Process 

The  MATLAB  codes  developed  in  this  research,  along  with  a  description  of  each 
code’s  function,  is  presented  in  detail  in  the  Appendix.  The  codes  themselves  contain 
comments  describing  the  purpose  of  every  section.  A  general  overview  of  the  software 
system  process  is  shown  here  in  Ligures  (13)  and  (14). 
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Figure  13.  Software  System  Process  Chart  Page  1  of  2 
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Figure  14.  Software  System  Process  Chart  Page  2  of  2 
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V.  Program  Validation 


Overview 

This  chapter  shows  comparisons  of  the  results  obtained  from  the  MATLAB 
reentry  program  with  those  obtained  from  other  established  software  programs  as  well  as 
the  results  that  are  predicted  by  first-order  analytical  approximations  of  the  equations  of 
motion.  First,  the  software  will  be  compared  with  two  trajectories  generated  by 
TARGET.  Next,  it  will  be  compared  to  a  trajectory  generated  by  IMPULSE.  The 
MATLAB  program  will  then  be  tested  against  the  first-order  solutions  at  near-circular 
orbital  velocities  for  a  shallow  gliding  entry,  a  steep  gliding  entry,  and  a  steep  ballistic 
entry.  Finally,  the  program  will  be  examined  against  the  first-order  approximation  for 
medium  gliding  entry  at  supercircular  orbital  velocity. 

Comparison  with  TARGET 

As  mentioned  in  the  background  for  this  research,  TARGET  was  developed  as  a 
tool  for  initial  analysis  of  entry  trajectories.  It  was  intended  to  be  used  as  a  substitute  for 
POST  during  the  systems  concept  stage  when  simplicity  and  speed  is  preferred  and  a 
good  estimate  will  suffice.  As  such,  TARGET  was  validated  against  POST  and  the 
trajectories  they  produced  agreed  quite  nicely.  TARGET  results  only  differed  from 
POST  by  less  than  .004%  for  flight-path  angle,  .001%  for  latitude  and  longitude,  and  .6% 
for  velocity  (Hicks,  1993:4-7).  By  validating  the  MATLAB  program  results  against 
TARGET,  it  is  also  implicitly  validated  against  POST. 
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Case  #1. 


The  first  comparison  consists  of  the  trajectory  for  a  notional  high-speed,  ballistic 
projectile.  Both  TARGET  and  the  MATLAB  program  begin  with  the  same  initial 
parameters  as  shown  in  Table  ( 1 ).  For  this  case  the  classic  orbital  elements  of  the  reentry 
vehicle  are  used  as  initial  conditions  so  the  function  of  the  MATLAB  program  to  convert 
these  elements  as  described  in  Chapter  III  is  tested.  TARGET  assumes  the  Earth’s 
atmosphere  is  spherical,  therefore  the  MATLAB  program  used  a  spherical  atmosphere  to 
remain  consistent.  The  trajectories  generated  by  each  program  are  plotted  in  Figures  (10) 
-(15). 


Table  (1):  Parameters  for  Comparison  with  TARGET  case  #1 


Mass:  13  kg 
CL:  0.00 

CD :  .08575 

S:  .036609615  m2 

flf(0):  -1840738  m 

v(0) :  4.552029  rad 
z(0) :  1.178101 

e(0):  1.889431 

co( 0):  2.377746  rad 
Q(0) :  4.213107  rad 
cr :  0  rad 

Corresponding  initial  conditions  for  equations 
of  motion  (automatically  generated) 

r(0):  6774400  m 

0(0):  4.4943  rad 

^(0):  .5903  rad 

RV (0) :  18220  m/s 

y( 0):  -1.2217  rad 

^•(0) :  1.1506  rad 
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West  Longitude  (deg) 


Figure  15.  Altitude/Time  Comparison  for  MATLAB/TARGET  Case  #1 


Altitude  (km) 


Figure  16.  Longitude/Altitude  Comparison  for  MATLAB/TARGET  Case  #1 
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Relative  Velocity  (km/sec)  North  Latitude  (deg) 


Figure  17.  Latitude/ Altitude  Comparison  for  MATLAB/TARGET  Case  #1 


Figure  18.  Velocity /Altitude  Comparison  for  MATLAB/TARGET  Case  #1 
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Heading  Angle  (deg)  Relative  Flight-Path  Angle  (deg) 


Altitude  (km) 

igure  19.  Flight-Path  Angle/ Altitude  Comparison  for  MATLAB/TARGET  Case  #1 


Figure  20.  Heading/ Altitude  Comparison  for  MATLAB/TARGET  Case  #1 
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In  order  to  quantify  the  results  of  the  comparison,  the  differences  were  summarized  on  a 
percentage  basis  and  plotted  against  time  in  Figures  (16)  and  (17).  The  percentage  values 
were  computed  at  each  interval  of  time  using: 


%  difference  = 


(TARGET  Result)  -  (MATLAB  Result) 
MATLAB  Result 


Figure  21.  Comparison  of  Angular  Data  for  MATLAB/TARGET  Case  #1 


(68) 
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7 


Figure  22.  Comparison  of  Velocities  for  MATLAB/TARGET  Case  #1 


These  results  illustrate  that  the  MATLAB  program  and  TARGET  produce 
comparable  results  when  starting  with  the  same  initial  parameters.  It  also  demonstrates 
the  validity  of  the  routine  used  in  MATLAB  to  convert  the  orbital  elements  to  the 
parameters  used  in  the  equations  of  motion.  The  deviation  in  velocity  at  the  lower  levels 
of  the  atmosphere  is  almost  certainly  a  result  of  the  different  atmospheric  models  used  by 
the  two  programs.  TARGET  uses  a  more  complex  seven-layered  atmosphere  while  the 
MATLAB  program  uses  a  simplified  single  layer  exponential  atmosphere. 


Case  #2. 

A  second  comparison  was  perfonned  between  the  MATLAB  routine  and 
TARGET  using  the  vehicle  tested  by  Albrecht  in  her  analysis  of  the  effects  of  lift  and 
thrust  on  maneuvering  a  reentry  vehicle  (2005:24-29).  This  vehicle  is  a  much  larger  and 
heavier  model  than  the  one  in  the  previous  test  and  includes  a  substantial  amount  of  lift  at 
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a  significant  bank  angle.  It  also  begins  its  trajectory  at  a  much  lower  velocity.  The  initial 
parameters  are  listed  in  Table  (2)  and  the  results  are  shown  in  Figures  (18)  -  (26). 


Table  (2):  Parameters  for  Comparison  with  TARGET  case  #2 


Mass: 

700  kg 

CL: 

.2 

CD : 

.2 

S: 

,216419m2 

*0): 

6500140  m 

0(0): 

-1.34732  rad 

m- 

.65176  rad 

*F(0):  3007.2  m/s 
y(0) :  -.523599  rad 
y/( 0):  1.195427  rad 
cr :  .8  rad 


Figure  23.  Altitude/Time  comparison  for  MATLAB/TARGET  Case  #2 
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North  Latitude  (deg)  West  Longitude  (deg) 
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Figure  24.  Longitude/Altitude  Comparison  for  MATLAB/TARGET  Case  #2 


Figure  25.  Latitude/ Altitude  Comparison  for  MATLAB/TARGET  Case  #2 
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Relative  Flight-Path  Angle  (deg)  Relative  Velocity  (km/sec) 


Figure  26.  Velocity /Altitude  Comparison  for  MATLAB/TARGET  Case  #2 


Figure  27.  Flight-Path  Angle/ Altitude  Comparison  for  MATLAB/TARGET  Case  #2 
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Figure  28.  Heading/ Altitude  Comparison  for  MATLAB/TARGET  Case  #2 
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Figure  29.  Comparison  of  Angular  Data  for  MATLAB/TARGET  Case  #2 
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Percent  Difference  in  Velocity  Difference  in  Flight-Path  Angle  (deg) 


Figure  30.  Comparison  of  Flight-Path  Angles  for  MATLAB/TARGET  Case  #2 


Figure  31.  Comparison  of  Velocities  for  MATLAB/TARGET  Case  #2 
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These  results  demonstrate  the  ability  of  the  MATLAB  routine  to  competently 
simulate  modeled  trajectories  for  vehicles  with  lift  and  a  bank  angle.  Once  again  the 
results  for  predicted  position  from  the  MATLAB  program  and  TARGET  are  virtually 
identical.  While  not  excessive,  there  are  deviations  in  the  predicted  velocity,  flight-path 
angle  and  heading  angle  between  the  two  programs.  This  is  probably  attributable  to  the 
different  atmospheric  models  used.  The  fact  that  Case  #2  was  a  high-lift  model  while 
Case  #1  was  a  ballistic  model  is  almost  certainly  the  reason  those  deviations  were  more 
pronounced  in  this  simulation  than  in  the  previous  one. 

Comparison  with  IMPULSE 

The  MATLAB  program  was  also  compared  to  IMPULSE  to  determine  if  they 
produced  similar  results.  As  described  in  the  Introduction,  IMPULSE  is  designed  to 
simulate  the  trajectories  of  ballistic  missiles  from  launch  to  impact.  Since  the  MATLAB 
program  does  not  consider  thrust,  the  comparison  is  made  starting  from  the  point 
immediately  following  burnout  and  separation  of  the  reentry  vehicle. 

One  problem  in  developing  the  comparison  between  the  two  programs  is  that  the 
output  from  IMPULSE  displays  the  latitude  in  geodetic  coordinates  rather  the  geocentric 
coordinates  used  in  Equations  (9)  -  (11)  and  (29)  -  (31).  Therefore,  the  geocentric 
latitude  output  from  the  MATLAB  program  was  converted  to  geodetic  latitude  for 
comparison  purposes  (Geocentric  2005). 

To  begin  with,  the  angle  shown  as  (f>*a  in  Figure  (27)  is  the  relationship  between 

the  geocentric  latitude  and  the  geodetic  latitude  at  the  point  where  r  crosses  the  Earth’s 
surface,  described  by 
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.*  tan  (b 

tan  t  =7, — 72 
(1-e) 


(69) 


where  e  is  the  ellipticity  of  the  Earth  described  in  Equation  (41).  The  angular  difference 
between  geocentric  and  geodetic  latitude  at  this  point  on  the  Earth  is: 

#  =  (70) 

The  radius  of  curvature  of  the  Earth  ( pa )  at  this  point  in  the  Meridian  is: 


Figure  32.  Reference  Ellipsoid 
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The  difference  between  (fa  and  the  geodetic  latitude  of  the  reentry  vehicle  ( </>*)  is  then 

approximated  as  a  function  of  the  height  above  the  surface  of  the  planet  where  h  and  /  in 
Figure  (5)  are  assumed  to  be  equal  for  the  Earth  as  expressed  in  Equation  (43): 


tan  8<j>  = 


h  sin  S(j) 
Pa+k 


Subtracting  8<f  from  <f>*a  gives  the  geodetic  latitude: 


*  r»  ,* 


4  =  &-<¥ 


(72) 


(73) 


The  vehicle  parameters  and  the  conditions  following  burnout,  described  in  the 
ECI  frame,  are  shown  in  Table  (3).  The  results  of  the  comparison  are  shown  in  Figures 
(28) -(35). 


Table  (3):  Parameters  for  Comparison  with  IMPULSE 


Mass: 

1088.64  kg 

Corresponding  initial  conditions  for  equations 
of  motion  (automatically  generated) 

CL: 

0.00 

CD : 

.11 

r(0): 

6946000  m 

S: 

1.880022  m2 

m-- 

-1.374  rad 

X(0): 

1044303.78549  m 

m-- 

.772  rad 

Y(0): 

-4866124.54174  m 

V(0): 

5194.5  m/s 

Z(0): 

4845298.63761  m 

r{  o): 

1.0775  rad 

X(0): 

-1094.1974  m/s 

K°): 

-2.7926 

Y(0): 

-4187.85837  m/s 

Z(0): 

2588.77299  m/s 

cr : 

0  rad 

t(0): 

00:03:20.6  (sidereal  time) 
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West  Longitude  (deg)  Altitude  (km) 


Figure  33.  Altitude/Time  Comparison  for  MATLAB/IMPULSE 


Figure  34.  Longitude/ Altitude  Comparison  for  MATLAB/IMPULSE 
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Velocity  (km/s)  North  Latitude  (deg) 


Figure  35.  Geodetic  Latitude/ Altitude  Comparison  for  MATLAB/IMPULSE 


Altitude  (km) 


Figure  36.  Velocity/Altitude  Comparison  for  MATLAB/IMPULSE 
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Figure  37.  Flight-Path  Angle/ Altitude  Comparison  for  MATLAB/IMPULSE 


Figure  38.  Comparison  of  Geographical  Angular  Data  for  MATLAB/IMPULSE 
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Percent  Difference  Difference  (deg) 


Figure  39.  Comparison  of  Flight-Path  Angles  for  MATLAB/IMPULSE 
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Figure  40.  Comparison  of  Velocities  for  MATLAB/IMPULSE 


55 


The  angular  data  results  are  essentially  the  same  for  both  programs.  The  fact  that 
the  initial  conditions  matched  validates  the  routine  used  in  MATLAB  to  convert  inertial 
coordinates  to  the  parameters  used  in  the  equations  of  motion.  The  difference  in  flight- 
path  angles  was  not  plotted  as  a  percent  difference  due  to  the  discontinuity  as  the  angles 
approached  zero.  The  flight-path  angle  and  altitude  comparisons  demonstrate  the  ability 
of  the  MATLAB  program  to  simulate  an  orbital  trajectory  far  from  the  planet’s  surface. 
Although  the  velocities  match  throughout  the  exoatmospheric  portion  of  the  trajectories, 
there  is  a  significant  deviation  as  the  reentry  vehicle  penetrates  the  denser  levels  of  the 
atmosphere.  This  is  probably  a  result  of  the  fact  that  IMPULSE  varies  CD  and  CL  as  a 
function  of  velocity  and  small  angles  of  pitch  and  yaw  while  the  MATLAB  program 
assumes  CD  is  constant  and  used  CL=  0  for  this  simulation.  Further  analyis  of  this 
deviation  in  follow-on  research  may  be  appropriate. 

Comparison  with  First-Order  Analytic  Solutions 

The  equations  of  motione  described  in  Chapter  I  can  be  simplified  and  solved 
analytically  for  special  cases  when  certain  assumptions  are  made.  Such  solutions  are 
approximations  but  provide  a  good  estimate  as  long  as  the  parameters  fit  the  assumptions 
made  for  each  special  case. 

Shallow,  Gliding  Entry. 

The  first  analytic  solution  tested  is  the  analysis  for  shallow,  gliding  entry.  This 
solution  assumes  the  vehicle  produces  enough  lift  to  maintain  hypersonic  glide  at  a  small 
flight-path  angle.  Once  the  flight-path  angle  exceeds  the  point  where  the  small  angle 
approximations 
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sin  y  x  y 

(74) 

cosy  ~  1 

(75) 

no  longer  hold  true,  the  analysis  is  no  longer  valid. 

The  solved  equations  of  motion  for  shallow,  gliding  entry  are  given  by  (Hicks, 
2005:81-86) 


pSCD 
2m  (3 


(76) 


where  r)  is  a  convenient  way  to  express  the  altitude  non-dimensionally  as  a  function  of 
density.  The  relative  velocity  is  related  to  the  altitude  by: 


Rrr2 


V  =  ■ 


S  oro 


fc  A 

\CdJ 


+  1 


(77) 


The  flight-path  angle  is  given  by: 


y, 


-2go 


fC  A 


1 3RV 2 


Ignoring  the  gravity  force  along  the  flight  path,  the  acceleration  is: 


dRV 

dt 


=  -/3ilRV 2 


(78) 


(79) 
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The  input  parameters  for  the  comparison  with  the  shallow,  gliding  entry  analytic 
solution  are  intended  to  simulate  a  large  high-drag  vehicle  with  moderate  lift  and  are 
given  in  Table  (4).  The  graphical  comparisons  are  shown  in  Figures  (36)  -  (40). 


Table  (4):  Parameters  for  Shallow,  Gliding  Entry  Comparison 


Mass: 

99000  kg 

CL: 

.5 

CD : 

1.0 

S: 

20  m2 

r(0): 

6438140  m 

r(  0): 

-.01304  rad 

RV  (0): 

7571.7  m/s 

m- 

.5236  rad 

m- 

0  rad 

HO): 

.64649  rad 

<j : 

0  rad 

Figure  41.  Velocity/Altitude  Comparison  for  Shallow,  Gliding  Entry 
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Figure  42.  Flight-Path  Angle/ Altitude  Comparison  for  Shallow,  Gliding  Entry 


Figure  43.  Deceleration/ Altitude  Comparison  for  Shallow,  Gliding  Entry 
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Figure  44.  Comparison  of  Velocities  for  Shallow,  Gliding  Entry 


Figure  45.  Comparison  of  Flight-Path  Angles  for  Shallow,  Gliding  Entry 
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These  results  validate  that  the  MATLAB  program  has  been  coded  correctly  as 
they  match  very  well  in  the  portion  of  the  trajectory  where  the  first-order  solution  is 
accurate.  Not  surprisingly,  the  first-order  solution  fails  to  properly  predict  the  vehicle’s 
behavior  when  the  small-angle  approximations  for  y  are  no  longer  valid,  as  illustrated 
particularly  well  in  Figure  (37).  The  predicted  velocity,  flight-path  angle,  and 
deceleration  all  begin  to  diverge  when  y  decreases  below  about  -3° ,  as  shown  in  Figures 
(38) -(40). 

Since  the  deviations  in  the  comparison  with  first-order  theory  were  obviously  the 
result  of  the  breakdown  of  the  small-angle  approximations  for  y,  a  quick  look  at  this 
problem  using  a  second-order  solution  is  appropriate.  Loh’s  unified  theory  for  entry 
trajectories  reveals  that  the  relative  flight-path  angle  is  given  by  (Hicks,  2006:148): 


C 


cosye+  '  {v-Ve) 


cos  y  =  - 


C 


1  + 


Ac 


o  V 


8oro 

Ry2 


-1 


1- 


% 

V  ) 


(80) 


Figures  (41)  and  (42)  show  the  results,  once  again  using  the  parameters  shown  in 
Table  (4).  The  flight-path  angle  predicted  by  the  second-order  approximations  matches 
the  angle  predicted  by  the  MATLAB  program  much  more  closely  than  the  first-order 
solution,  especially  as  the  flight-path  angle  decreases  below  -3°.  These  results  further 
bolster  confidence  in  the  MATLAB  routine  and  help  explain  the  deviations  encountered 
with  first-order  theory  as  the  flight-path  angle  steepens. 
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Figure  46.  Second-Order  Flight-Path  Angle  Comparison  for  Shallow,  Gliding  Entry 


Figure  47.  Second-Order  Comparison  of  Flight-Path  Angles  for  Shallow,  Gliding  Entry 
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Steep,  Gliding  Entry. 

The  second  comparison  with  analytical  solutions  tests  the  case  of  steep  entry  with 
lift  at  near-circular  velocity.  Using  the  assumption  that  drag  and  lift  are  the  dominant 
factors  affecting  the  trajectory,  Hicks  derives  the  approximate  solutions  for  relative 
velocity  and  relative  flight-path  angle  and  deceleration  ( adecel )  as  (2005:55-57): 
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The  subscript  e  denotes  the  initial  condition  at  entry. 

The  model  used  for  this  test  is  a  small,  dense,  streamlined  vehicle  with  lift 
entering  at  a  near-vertical  flight-path  angle.  The  initial  parameters  are  listed  in  Table  (5) 
and  the  results  are  shown  in  Figures  (43)  -  (47). 
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Table  (5):  Parameters  for  Steep,  Gliding  Entry  Comparison 


Mass: 

50  kg 

CL: 

.03 

CD : 

.01 

S: 

.01  m2 

r(0): 

6478140  m 

r(  0): 

-85° 

RV  (0): 

5372  m/s 

m- 

.600588  rad 

0(0): 

-1.8075  rad 

K0): 

1.916  rad 

cr : 

0  rad 

Figure  48.  Velocity/Altitude  Comparison  for  Steep,  Gliding  Entry 
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Flight-Path  Angle  (deg) 


Figure  49.  Flight-Path  Angle/ Altitude  Comparison  for  Steep,  Gliding  Entry 


Figure  50.  Deceleration/ Altitude  Comparison  for  Steep,  Gliding  Entry 
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Figure  51.  Comparison  of  Data  for  Steep,  Gliding  Entry 


Figure  52.  Comparison  of  Deceleration  difference  for  Steep,  Gliding  Entry 
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The  close  agreement  of  these  results  strongly  suggest  that  the  equations  have  been 
coded  correctly  in  MATLAB  and  that  the  program  produces  reliable  trajectories  for 
steep,  gliding  entry  problems. 

Steep,  Ballistic  Entry. 

First-Order  approximation  for  steep,  ballistic  entry  models  a  hypersonic  vehicle 
with  no  lift  at  a  near-vertical  flight-path  angle.  The  solution  makes  the  assumption  that 
the  flight-path  angle  remains  constant  throughout  the  trajectory.  The  simplified  equations 
of  motion  for  relative  velocity  and  relative  flight-path  angle  are  derived  by  Hicks 
(2006:118-121). 

Since  the  flight-path  angle  is  a  constant: 


Y  =  Ye 


(84) 


For  simplification  the  variable  a  is  defined  as 


a  = 


-2/7 
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and  a  constant  A  defined  as 
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where  Ei ( ae )  is  the  exponential  integral  ae.  The  relative  velocity  is  defined  as: 
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Once  again,  the  subscript  e  denotes  the  initial  entry  condition. 

The  initial  parameters  for  the  comparison  between  the  MATLAB  program  and 
first-order  analytical  solution  for  steep,  ballistic  entry  are  shown  in  Table  (6)  and  the 
results  are  shown  graphically  in  Figures  (48)  -  (50). 


Table  (6):  Parameters  for  Steep,  Ballistic  Entry  Comparison 


Mass: 

13  kg 

CL: 

0 

CD : 

.08575 

S: 

.036609615  m2 

r(0): 

6778140  m 

7(0): 

-85° 

RV  (0): 

10193  m/s 

m- 

.58059  rad 

0(0): 

-1.8875  rad 

V(  0): 

1.7  rad 

<j : 

0  rad 
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Figure  53.  Velocity /Altitude  Comparison  for  Steep,  Ballistic  Entry 


Figure  54.  Flight-Path  Angle/ Altitude  Comparison  for  Steep,  Ballistic  Entry 
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Figure  55.  Comparison  of  Data  for  Steep,  Ballistic  Entry 


Once  again,  the  first-order  theory  matches  well  with  the  MATLAB  program.  The 
small  deviation  in  flight-path  angle  is  expected  because  the  analytical  solution  treats  this 
parameter  as  a  constant,  as  expressed  in  Equation  (84).  The  velocity  difference  is 
consistent  with  the  previous  comparisons  with  first-order  solutions. 


Medium,  Gliding  Entry  at  Supercircular  Velocity. 

The  final  test  is  another  look  at  medium,  gliding  entry,  but  this  time  there  is  no 
restriction  that  the  initial  velocity  must  be  at  near-circular  speed.  Hicks  shows  how 
Wang  and  Ting  developed  their  approximate  solution  for  this  condition  with  the 
assumption  that  the  flight-path  angle  remains  relatively  small  (2006:104-107).  Their 
simplified  equations  for  flight-path  angle  and  relative  velocity  are  as  follows: 
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As  with  shallow,  gliding  entry  at  near-circular  speeds,  the  acceleration  remains 


dRV 

dt 


=  -prjRV 1 


(90) 


The  simulated  model  for  the  comparison  with  the  first-order  solution  for  medium, 
gliding  entry  at  supercircular  velocity  is  a  small  lifting  body  entering  the  atmosphere  at  a 
moderate  flight-path  angle  with  high  velocity.  The  initial  parameters  are  listed  in  Table 
(7)  and  the  results  are  shown  in  Figures  (5 1)  -  (54). 


Table  (7):  Parameters  for  Medium,  Gliding  Entry  at  Supercircular  Velocity 
Comparison 


Mass: 

2000  kg 

CL: 

.075 

CD : 

.3 

S: 

20  m2 

r(0): 

6478140  m 

r(  0): 

-35° 

RV  (0): 

9847  m/s 

m- 

.5236  rad 

0(0): 

0  rad 

AO): 

1.146  rad 

<j : 

0  rad 
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Figure  56.  Velocity /Altitude  Comparison  for  Medium,  Gliding  Entry  at  Supercircular  Speed 
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Figure  57.  Flight-Path  Angle/Altitude  Comparison  for  Medium,  Gliding  Entry  at  Supercircular 

Speed 
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Figure  58.  Deceleration/ Altitude  Comparison  for  Medium,  Gliding  Entry  at  Supercircular  Speed 


Figure  59.  Comparison  of  Data  for  Medium,  Gliding  Entry  at  Supercircular  Speed 
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Due  to  the  fact  that  the  first-order  solution  for  flight-path  angle  becomes 
imaginary  when  the  quantity  in  brackets  on  the  right  side  of  Equation  (88)  goes  below 
zero,  the  comparison  could  only  be  plotted  to  an  altitude  of  2 1 .3  km.  Below  this  altitude 
the  first-order  solutions  for  flight-path  angle,  relative  velocity,  and  acceleration  became 
complex  numbers  and  were  no  longer  valid  for  comparison  purposes.  Up  to  the  point 
where  the  analytic  solution  fell  apart,  the  MATLAB  program  agreed  closely  with  the 
first-order  approximations,  further  solidifying  the  fact  that  the  program  is  operating 
properly. 
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VI.  Conclusion  and  Recommendations 


Conclusion 

For  this  research  an  algorithm  was  developed  using  MATLAB  to  model  a 
trajectory  for  a  vehicle  entering  the  Earth’s  atmosphere.  The  algorithm  used  the 
equations  of  motion  developed  in  Chapter  II  and  established  the  initial  conditions  using 
the  methods  shown  in  Chapter  III  to  establish  the  simulation.  The  program  was  then 
compared  to  previously  established  software  that  were  also  designed  to  predict  entry 
trajectories.  It  was  also  compared  to  first-order  approximations  for  specific  types  of 
reentry.  The  program  code  is  listed  in  the  Appendix. 

The  MATLAB  code  generated  fulfills  the  goal  of  this  research  of  creating  a 
program  to  quickly  generate  atmospheric  entry  simulations.  It  compares  very  favorably 
with  TARGET,  a  FORTRAN  program  with  similar  objectives,  and  IMPULSE,  another 
MATLAB  program  focused  on  simulating  ICBMs.  In  both  cases,  the  methods  used  in 
Chapter  III  to  establish  the  initial  conditions  yielded  perfect  results,  confirming  that  the 
equations  used  were  valid  and  coded  correctly. 

Comparison  with  the  first-order  approximations  further  solidified  the  validity  of 
the  MATLAB  algorithm.  The  MATLAB  program  routinely  produced  results  that  closely 
matched  the  first-order  solutions  over  those  portions  of  the  trajectory  one  would  expect 
the  analytic  approximations  to  maintain  accuracy. 

Although  the  model  used  in  this  research  is  not  all  encompassing,  it  provides  the 
capability  to  produce  reliable  results  for  systems-level  studies.  The  program  was 
intended  to  be  simple  and  quick  for  initial  analysis  and  to  utilize  MATLAB  due  to  the 
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language’s  emerging  widespread  employment.  The  algorithm  is  not  intended  to  be  a 
replacement  for  more  comprehensive  routines  that  already  exist,  such  as  POST. 

Recommendations  for  Future  Work 

The  MATLAB  routine  developed  in  this  research  lends  itself  to  future 
development.  There  are  several  areas  of  the  algorithm  outside  the  scope  of  this  research 
that  could  be  improved.  The  software  is  quite  adaptable  and  implementing  future 
adaptations  should  not  be  too  difficult. 

This  research  used  a  simple  exponential  model  as  an  adequate  approximation  of 
the  Earth’s  atmosphere.  Improvement  could  be  made  by  more  completely  modeling  the 
atmosphere  as  several  layers,  each  with  a  different  scale  height. 

For  simplification,  this  research  made  the  assumption  that  CL  and  CD  are 

constant  throughout  the  atmospheric  entry.  In  reality,  both  coefficients  are  functions  of 
angle  of  attack  and  velocity.  An  algorithm  could  be  implemented  to  allow  changes  in 
angle  of  attack  at  given  altitudes  or  velocities  and  to  calculate  the  corresponding 
coefficients  of  lift  and  drag. 

Just  as  angle  of  attack  was  held  constant,  the  bank  angle  (cr)  was  also  assumed  to 

be  constant  during  entry.  By  allowing  variations  in  cr  and  angle  of  attack,  a  future 
development  of  the  software  would  be  able  to  simulate  various  trajectories  that  use  lift  to 
maneuver  the  vehicle  during  reentry. 

This  research  also  assumed  that  the  reentry  vehicle  had  no  propulsion  system. 
The  equations  of  motion  as  shown  by  Hicks  (2006:58)  account  for  a  thrust  component. 
These  equations  could  be  substituted  for  the  ones  shown  in  Chapter  II  and  adapted  to 
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simulate  a  vehicle  that  uses  thrust  or  a  combination  of  thrust  and  lift  to  maneuver  during 
reentry. 

Finally,  future  research  may  be  interested  in  the  amount  of  aerodynamic  heating 
expected  during  reentry.  Vinh  (1980:141-156)  shows  how  heating  rates  vary  with 
velocity  and  atmospheric  density.  These  analytic  solutions  could  be  added  to  the  existing 
code  to  generate  the  heating  information  for  the  entire  trajectory. 

In  all  of  these  cases,  the  MATLAB  program  generated  in  this  research  would  be 
the  core  program.  The  routine  would  maintain  the  ability  to  generate  quick  and  simple 
results  with  a  reasonably  high  degree  of  accuracy  for  high-level  studies. 


77 


Appendix:  MATLAB  Code 


Core  Code 

This  is  the  MATLAB  core  code  for  simulating  entry  trajectories  (reentry. m).  It 
uses  the  algorithms  developed  in  this  research.  The  codes  for  operating  the  graphical 
user  interfaces  (inputs. m,  vpoint.m,  eci.m  and  orb.rri)  and  the  ordinary  differential 
equation  solvers  (spherical. m,  ellipsoid.m )  are  listed  following  the  core  code. 


%  Capt  Robert  Jameson 
%  AFIT/GSS/ENY/06-08 
%  14  January  2006 

%  This  program  simulates  the  reentry  trajectory  of  a  vehicle  with  perfect 
%  knowledge  of  the  vehicle  parameters  and  the  initial  conditions  for  the 
%  equations  of  motion. 


clc 

%  Clear  acceleration  data 
clear  'A' 

%  Surface  density  of  atmosphere  in  kg/mA3  at  STP  (assumption) 
rhos=l. 54793; 

%  Average  scale  height  in  meters  (assumption) 

H=6935; 

%  Calculate  inverse  of  scale  height  (beta) 
beta=l/H; 

%  Radius  of  the  Earth  at  the  equator  in  meters  (WGS-84  value) 
re=6378137; 

%  Radius  of  the  Earth  at  the  poles  in  meters  (WGS-84  value) 
rp=6356752.3142; 

%  This  equation  calculates  the  eccentricity  of  the  Earth's  atmosphere  (ee) 
ee=l-rp/re; 

%  WGS-84  acceleration  due  to  gravity  at  the  radius  of  equator  (m/sA2) 
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g0=9.798; 


%  WGS-84  rotation  rate  of  Earth  (rad/s) 
omegae=. 000072921 151467; 

%  Create  initial  user  inputs  if  none  exist 
if  ~exist(’user_inputs') 

%  Run  GUI  (inputs. m)  to  get  user  inputs  for  the  shape  of  the  Earth's 
%  atmosphere  and  the  desired  coordinate  system 
inputs 

uiwait(inputs) 

%  Save  values  for  coordinate  system  (l=vehicle  pointing,  2=ECI, 

%  3=classic  orbital  elements)  and  shape  of  atmosphere  (l=spherical, 
%  2=ellipsoid)  into  .mat  files 
load(’var') 
load('shape') 

%  If  the  vehicle  pointing  system  is  selected,  run  GUI  for  user  inputs  for 
%  vehicle  parameters  and  initial  conditions  (vpoint.m) 
if  var==  1 
vpoint 

uiwait(vpoint) 

%  Save  values  for  vehicle  parameters  (m=mass,  Cl=coefficient  of  lift, 
%  Cd=coefficient  of  drag,  S=reference  surface  area)  and  initial 
%  conditions  (ri=radius,  gammar=relative  flight-path  angle, 

%  vr=relative  velocity,  thetai=longitude,  phii=geocentric  latitude, 

%  psir=relative  heading,  sigma=bank  angle)  into  .mat  files 
load(’m') 
load(’Cl') 
load('Cd’) 
load(’S') 
load(’ri’) 
load(’gammar’) 
load('vr') 
load(’thetai’) 
load(’phii’) 
load('psir') 
load(’sigma') 
end 

%  If  the  ECI  system  is  selected,  run  GUI  for  user  inputs  for 
%  vehicle  parameters  and  initial  conditions  (eci.m) 
if  var==2 
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eci 

uiwait(eci) 

%  Save  values  for  vehicle  parameters  (m=mass,  Cl=coefficient  of  lift, 

%  Cd=coefficient  of  drag,  S=reference  surface  area)  and  initial 
%  conditions  (x,y,z=inertial  position,  vx,vy,vz=inerial  velocity, 

%  sigma=bank  angle,  hour, min, sec=Sidereal  time)  into  .mat  files 
load('m') 
load('Cl') 
load('Cd') 
load(’S') 
load(’x’) 
load(’y') 
load(’z') 
load(’vx’) 
load('vy') 
load('vz’) 
load(’sigma’) 
load(’hour') 
load(’min') 
load(’sec') 
end 

%  If  the  classic  orbital  elements  are  selected,  run  GUI  for  user  inputs  for 
%  vehicle  parameters  and  initial  conditions  (orb.m) 
if  var==3 
orb 

uiwait(orb) 

%  Save  values  for  vehicle  parameters  (m=mass,  Cl=coefficient  of  lift, 

%  Cd=coefficient  of  drag,  S=reference  surface  area)  and  initial 
%  conditions  (ri=radius,  gammar=relative  flight-path  angle, 

%  vr=relative  velocity,  thetai=longitude,  phii=geocentric  latitude, 

%  psir=relative  heading,  sigma=bank  angle)  into  .mat  files 
load(’m') 
load(’Cl’) 
load('Cd') 
load(’S’) 
load(’ri') 
load('gammar’) 
load('vr') 
load(’thetai’) 
load('phii’) 
load('psir’) 
load(’sigma') 
end 
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%  If  the  ECI  system  is  selected,  run  GUI  for  user  inputs  for 
%  vehicle  parameters  and  initial  conditions  (eci.m) 
if  var==2 
eci 

uiwait(eci) 

%  Save  values  for  vehicle  parameters  (m=mass,  Cl=coefficient  of  lift, 

%  Cd=coefficient  of  drag,  S=reference  surface  area)  and  initial 
%  conditions  (x,y,z=inertial  position,  vx,vy,vz=inerial  velocity, 

%  sigma=bank  angle,  hour, min, sec=Sidereal  time)  into  .mat  files 
load('m') 
load('Cl’) 
load(’Cd') 
load(’S’) 
load(’x’) 
load(’y’) 
load('z’) 
load('vx') 
load(’vy’) 
load(’vz’) 
load('sigma') 
load('hour’) 
load('min') 
load(’sec') 
end 

%  If  the  classic  orbital  elements  are  selected,  run  GUI  for  user  inputs  for 
%  vehicle  parameters  and  initial  conditions  (orb.m) 
if  var==3 
orb 

uiwait(orb) 

%  Save  values  for  vehicle  parameters  (m=mass,  Cl=coefficient  of  lift, 

%  Cd=coefficient  of  drag,  S=reference  surface  area)  and  initial 
%  conditions  (a=semi-major  axis,  nu=true  anomaly,  inc=inclination, 

%  e=eccentricity,  omega=argument  of  perigee,  node=right  ascension  of 
%  ascending  node,  sigma=bank  angle,  hour, min, sec=Sidereal  time)  into  .mat 
%  fdes 

load('m') 

load('Cl’) 

load('Cd') 

load(’S’) 

load(’a’) 

load(’nu’) 

load('inc') 
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load(’e’) 

load('omega') 

load('node’) 

load(’sigma’) 

load(’hour') 

load('min') 

load(’sec') 

end 

%  set  up  userinputs  and  parameters  matrix  if  vehicle  pointing  system 
%  selected 
if  var==l 

user_inputs=[S;Cl;Cd;m;vr;garnmar;sigma;thetai;phii;ri;psir]; 

parameters=  {'Surface  Area';'Cl';'Cd';'Mass';'Entry  Velocity';... 
'Initial  Flight-Path  Angle’;... 

'Bank  Angle';'Longitude';’Latitude';’Initial  Radius';'Heading’}; 

%  Display  user  inputs  and  parameters  in  array  editor  for  user  reference 
%  and  modification 

openvar('parameters'); 

openvar('userinputs'); 

end 

%  set  up  user  inputs  and  parameters  matrix  if  ECI  system  selected 
if  var==2 

user_inputs=[S;Cl;Cd;m;x;y;z;vx;vy;vz;sigma;hour;min;sec]; 


parameters={ 'Surface  Area';'Cr;'Cd';'Mass’;'X’;’Y';... 

’Z’;’X  dot’;’Y  dot’;’Z  dof;’Bank  Angle';'Hour';’Minute';'Second’}; 

%  Display  user  inputs  and  parameters  in  array  editor  for  user  reference 
%  and  modification 

openvar('parameters'); 

openvar('userinputs’); 

end 

%  set  up  user  inputs  and  parameters  matrix  if  classic  orbital  elements 
%  selected 
if  var==3 

user_inputs=[S;Cl;Cd;m;a;e;inc;node;omega;nu;sigma;hour;min;sec]; 
parameters=  {'Surface  Area';'Cl';'Cd';'Mass';’Semi-Major  Axis';'Eccentricity';... 
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'Inclination';'Node';'Perigee';True  Aanomaly';'Bank 
Angle';'Hour';'Minute';’Second'}; 


%  Display  userinputs  and  parameters  in  array  editor  for  user  reference 
%  and  modification 

openvar('parameters’); 

openvar(’user_inputs'); 

end 

end 

%  Allow  editing  of  user  inputs  in  display  editor 
if  var==l 

S=user_inputs(l); 

Cl=user_inputs(2); 

Cd=user_inputs(3); 

m=user_inputs(4); 

vr=user_inputs(5); 

gammar=user_inputs(6); 

sigma=user_inputs(7); 

thetai=user_inputs(8); 

phii=user_inputs(9); 

ri=user_inputs(  1 0) ; 

psir=user_inputs(l  1); 

end 


if  var==2 


%  Allow  editing  of  user  inputs  in  display  editor 
S=user_inputs(l); 

Cl=user_inputs(2); 

Cd=user_inputs(3); 
m=user_inputs(4); 
x=user_inputs(5); 
y=user_inputs(6); 
z=user_inputs(7); 
vx=user_inputs(8); 
vy=user_inputs(9); 
vz=user_inputs(  1 0); 
sigma=user_inputs(  1 1); 
hour=user_inputs(  1 2); 
min=user_inputs(  13); 
sec=user_inputs(  1 4); 
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%  Develop  initial  conditions  when  ECI  system  is  selected 

%  Calculate  initial  radius  from  inertial  position  (ri) 
ri=(xA2+yA2+z  A2)A .  5 ; 

%  Calculate  magnitude  of  inertial  velocity  (vi) 
vi=(vxA2+vyA2+vzA2)  A .  5 ; 

%  Generate  vectors  for  inertial  position  and  velocity 
Rxyz=[x;y;z]; 

V  xyz=[vx;  vy ;  vz] ; 

%  Calculate  angular  momentum  vector  (H) 
H=cross(Rxyz,Vxyz); 

%  Calculate  magnitude  of  angular  momentum  (h) 
h=(H(  1  )A2+H(2)A2+H(3)A2)A.5; 

%  Calculate  the  inertial  flight-path  angle  (gammai)  and  find  the 

%  correct  sign. 

if  dot(Rxyz,Vxyz)>0 

gammai=ac  os(h/(ri  *  vi)) ; 
else 

gammai=-acos(h/(ri*vi)); 

end 

%  Calculate  the  geocentric  latitude  (phii) 
phii=asin(z/ri); 

%  Calculate  the  longitude  assuming  t=0  (theta) 
theta=asin(y/(xA2+yA2)A.5); 

%  Correct  longitude  for  sign  ambiguity 
if  x<0 

theta=pi-theta; 

end 

%  Correct  longitude  for  Greenwich  sidereal  time  (thetai) 
thetai=theta-(user_inputs(  1 2 , 1  )*3600+user_inputs(  1 3, 1  )*60+... 
user_inputs(  1 4, 1 ))  *  omegae; 

%  Ensure  the  magnitude  of  thetai  is  not  greater  than  2pi 
if  thetai<-2*pi 
thetai=thetai+2  *pi 
end 
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%  Calculate  the  inclination  of  the  orbit  (inc) 
inc=acos(H(3)/h); 

end 

%  Develop  initial  conditions  when  orbital  elements  are  selected 
if  var==3 


%  Allow  editing  of  user  inputs  in  display  editor 
S=user_inputs(l); 

Cl=user_inputs(2); 

Cd=user_inputs(3); 
m=user_inputs(4); 
a=user_inputs(5); 
e=user_inputs(6); 
inc=user_inputs(7); 
node=user_inputs(8); 
omega=user_inputs(9); 
nu=user_inputs(  1 0); 
sigma=user_inputs(  1 1); 
hour=user_inputs(  1 2); 
min=user_inputs(  13); 
sec=user_inputs(  1 4); 

%  Earth's  gravitational  parameter  (mA3/sA2) 
mu=3 .9860 1  *  1 0A 1 4; 

%  Calculate  the  semi-latus  rectum  of  the  orbit 
p=a*(l-eA2); 

%  Calculate  radius  (ri) 
ri=a*(  1  -eA2)/(  l+e*cos(nu)); 

%  Calculate  P-component  of  velocity  (vp)  in  the  perifocal  coordinate  system 
vp=(mu/p)A.5*(-sin(nu)); 

%  Calculate  Q-component  of  velocity  (vq)  in  perifocal  system 
vq=(mu/p)A.5*(e+cos(nu)); 

%  Establish  velocity  vector  in  perifocal  system 
V=[vp;vq]; 

%  Establish  radius  vector  in  perifocal  system 
R=[ri*cos(nu);ri*sin(nu)]; 
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%  Calculate  position  in  ECI  system  from  perifocal  system 

x=R(l)*(cos(node)*cos(omega)-sin(node)*sin(omega)*cos(inc))+R(2)*... 

(-cos(node)*sin(omega)-sin(node)*cos(omega)*cos(inc)); 

y=R(  1  )*  (sin(node)  *  cos(omega)+cos(node)  *  sin(omega)*  cos(inc))+R(2)  * . . . 
(-sin(node)  *  sin(omega)+cos(node)  *  cos(omega)*  cos(inc)); 

z=R(  1 )  *  (sin(omega)*  sin(inc))+R(2)  *  (cos(omega)  *  sin(inc)); 

%  Calculate  velocity  in  ECI  system  from  perifocal  system 

vx=V  ( 1  )*  (cos(node)*  cos(omega)-sin(node)  *  sin(omega)*  cos(inc))+ V  (2)  * . . 
(-cos(node)*sin(omega)-sin(node)*cos(omega)*cos(inc)); 

vy=V  ( 1  )*  (sin(node)  *  cos(omega)+cos(node)*  sin(omega)  *  cos(inc))+V  (2)  * . 
(-sin(node)*sin(omega)+cos(node)*cos(omega)*cos(inc)); 

vz=V  ( 1  )*  (sin(omega)  *  sin(inc))+V  (2)  *  (cos(omega)  *  sin(inc)); 

vi=(vxA2+vyA2+vzA2)  A .  5 ; 

%  Establish  radius  and  velocity  vectors  in  ECI  system 
Rxyz=[x;y;z]; 

Vxyz=[vx;vy;vz]; 

%  Calculate  angular  momentum  vector 
H=cross(Rxyz,Vxyz); 

%  Calculate  magnitude  of  angular  momentum 
h=(H(  1  )A2+H(2)A2+H(3)A2)A.5; 

%  Calculate  inertial  flight-path  angle  (gammai)  and  correct  for  sign 
%  ambiguity, 
if  nu>pi 

gammai=-acos(h/(ri*vi)); 

else 

gammai=ac  os(h/(ri  *  vi)) ; 
end 

%  Calculate  geocentric  latitude 
phi  i=asin(z/ri); 

%  Calculate  the  longitude  assuming  t=0  (theta) 
theta=asin(y/ (xA2 +yA2)  A .  5  ) ; 

%  Correct  for  ambiguity  of  theta 
if  x<0 
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theta=pi-theta; 

end 

%  Correct  longitude  for  Greenwich  sidereal  time  (thetai) 
thetai=theta-(user_inputs(  1 2 , 1  )*3600+user_inputs(  1 3 , 1 )  *  60+. . . 
user_inputs(  1 4, 1 ))  *  omegae; 

%  Ensure  the  magnitude  of  thetai  is  not  greater  than  2pi 
if  thetai<-2*pi 

thetai=thetai+2  *pi; 
end 

end 

%  Continue  determing  initial  conditions  when  ECI  system  or  classic  orbital 
%  elements  are  selected 
if  var>l 

%  Calculate  the  local  horizontal  component  (vh)  of  the  inertial  velocity 
vh=vi*cos(gammai); 

%  Calculate  the  local  velocity  of  the  Earth’s  rotation  (ve) 
ve=ri*omegae*cos(phii); 

%  Find  the  right  ascension  of  the  ascending  node  (node)  when  the  ECI 
%  system  is  selected 
if  var==2 

kH=cross([0;0;l],H); 

N=kH/((kH(l)A2+kH(2)A2+kH(3)A2)A.5); 

node=acos(N(l)); 

%  Correct  for  sign  ambiguity  of  right  ascension  of  ascending  node 
if  N(2)<0 
node=-node; 
end 


end 

%  Calculate  the  inertial  heading  (psii) 
psii=acos(cos(inc)/cos(phii)); 

%  Correct  for  sign  ambiguity  of  the  heading 
if  pi/2<abs(node-theta)  &  abs(node-theta)<3*pi/2 
psii=-psii; 

end 
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%  Calulate  the  value  c  from  law  of  cosines  (c  is  used  for  simplification 
%  in  the  calculation  of  relative  heading) 
c=(veA2+vhA2-2*ve*vh*cos(psii))A.5; 

%  Calculate  relative  heading  (psir) 
psir=acos((vh*cos(psii)-ve)/c); 

%  Correct  for  sign  ambiguity  of  relative  heading 
if  pi/2<abs(node-theta)  &  abs(node-theta)<3*pi/2 
psir=-psir; 
end 

%  Calculate  relative  velocity  (vr) 

vr=(viA2+veA2-2  *  vh*  ve*  cos(psii))A .  5 ; 

%  The  initial  flight  path  angle  relative  to  the  rotating  earth  (gammar) 

%  as  a  function  of  inertial  velocity,  inertial  flight  path  angle  and 
%  relative  velocity 

gammar=asin(  vi  *sin(  gammai)/  vr) ; 

end 

%  Establish  the  initial  condition  vector 
yO=[ri  vr  gammar  thetai  phii  psir]; 

%  Set  initial  and  final  time  for  integration.  Modify  tf  as  needed  to 
%  include  entire  trajectory 
ti=0; 

tf=  10000; 

%  Select  set  of  differential  equations  based  on  shape  of  atmosphere 
%  selected 
if  shape==l 

%  Options-yls.m  stops  integration  at  surface  of  spherical  Earth  (sea 
%  level),  initial  step  size  and  maximum  step  size  can  be  modified  as 
%  desired 

options=odeset('Events',@stops,’InitialStep',.0004,'MaxStep',l); 

[t  Y]=ode45(@spherical,[ti  tf],y0, [options], m,S, Cl, Cd, sigma); 

else 

%  Options-yl.m  stops  integration  at  surface  of  ellipsoid  Earth  (sea  level) 

options=odeset('Events’,@stope,'InitialStep',.0004,'MaxStep',l); 

[t  Y]=ode45(@ellipsoid,[ti  tf],y0, [options], m,S, Cl, Cd, sigma); 
end 
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%  Calculate  acceleration  of  vehicle  (A)  from  equation  of  motion  (dV) 
for  i=l:length(t) 

A(i,  1  )=-((rhos*exp(-beta*  (Y  (i,  1  )-re))*  S  *Cd*  Y  (i,2)A2)/2)/m-gO* . . . 
(re/Y(i,l))A2*sin(Y(i,3))+Y(i,l)*omegaeA2*cos(Y(i,3))*... 
(cos(Y  (i,5))*sin(Y  (i,3))-sin(Y  (i,5))*sin(Y  (i,6))*cos(Y  (i,3))); 
end 


if  shape==l 

%  Plot  Altitude/Time 

figure(  1  );plot(t,(  Y  ( : ,  1  )-re)/ 1000); 
xlabel('Time  (sec)’); 
ylabel(Altitude  (km)’); 

%  Plot  deceleration/Altitude 

figure(2);plot((Y(:,l)-re)/1000,-A(:,l)/g0); 
xlabel(’Altitude  (km)’); 
ylabel(’Deceleration/gO’); 

%  Plot  flight-path  angle/altitude 

figure(3);plot((Y(Cl)-re)/1000,Y(:,3)*180/pi); 

xlabel(’Altitude  (km)’); 

ylabel(’Relative  Flight-Path  Angle  (deg)’); 

%  Plot  Longitude/altitude 

figure(4);plot((Y(: ,  1  )-re)/l  000,  Y(:  ,4)*  1 80/pi); 
ylabel(’Longitude  (deg)’); 
xlabel(’Altitude  (km)’); 

%  Plot  Latitude/altitude 

figure(5);plot((Y(: ,  1  )-re)/ 1 000,  Y  (:  ,5)*  1 80/pi); 
ylabel(’North  Latitude  (deg)’); 
xlabel(’Altitude  (km)’); 

%  Plot  relative  velocity/altitude 

figure(6);plot((Y(: ,  1  )-re)/ 1 000,  Y(:, 2)/ 1000); 
xlabel(’Altitude  (km)’); 
ylabel(’Relative  Velocity  (km/sec)’); 

%  Plot  relative  heading/altitude 

figure(7);plot((Y(: ,  1  )-re)/ 1 000,  Y  (:  ,6)*  1 80/pi); 
ylabel(’Heading  Angle  (deg)’); 
xlabel(’Altitude  (km)’); 

else 

%  Plot  Altitude/Time 
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figure(  1  );plot(t,(Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000); 
xlabel('Time  (sec)'); 
ylabel(’ Altitude  (km)'); 

%  Plot  deceleration/Altitude 

figure(2);plot((Y(:,  l)-rp/(cos(Y  (length(t),5))* .. . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000,  -  A(: ,  1  )/g0); 
xlabel('Altitude  (km)’); 
ylabel(’Deceleration/gO'); 

%  Plot  flight-path  angle/altitude 

figure(3);plot((Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000, Y(:, 3)*  1 80/pi); 
xlabel(’Altitude  (km)'); 
ylabel(’Relative  Flight-Path  Angle  (deg)’); 

%  Plot  Longitude/altitude 

figure(4);plot((Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000,Y(:  ,4)*  1 80/pi); 
ylabel(’Longitude  (deg)'); 
xlabel(’Altitude  (km)'); 

%  Plot  Latitude/altitude 

figure(5);plot((Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y(length(t),5)))A2+(l-ee)A2)A.5))/1000,Y(:, 5)*  180/pi); 
ylabel('North  Latitude  (deg)'); 
xlabel('Altitude  (km)’); 

%  Plot  relative  velocity/altitude 

figure(6);plot((Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000,Y(:  ,2)/ 1 000); 
xlabel(’Altitude  (km)'); 
ylabel(’Relative  Velocity  (km/sec)'); 

%  Plot  relative  heading/altitude 

figure(7);plot((Y(: ,  1  )-rp/(cos(Y  (length(t),5))* . . . 

((tan(Y  (length(t),5)))A2+(  1  -ee)A2)A.5))/ 1 000, Y(:, 6)*  1 80/pi); 
ylabel(’Heading  Angle  (deg)'); 
xlabel('Altitude  (km)’); 

end 
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GUI  Files 


The  graphical  user  interface  (GUI)  files  allow  intuitive  and  simple  operation  of 
the  MATLAB  program.  They  are  called  from  the  core  program  as  needed. 

Coordinate  System  and  Atmosphere 

The  first  GUI  used  by  the  software  is  the  inputs. m  file.  This  file  allows  the  user  to 
select  a  coordinate  system  for  the  initial  known  parameters  from  a  list  including  the 
vehicle  pointing  system,  the  ECI  system  and  classic  orbital  elements.  It  also  allows  the 
user  to  specify  a  spherical  atmosphere  or  an  ellipsoidal  atmosphere.  The  GUI  display  is 
shown  in  Figure  (52)  followed  by  the  MATLAB  code  for  inputs. m. 


Figure  60:  inputs.m  GUI  Display 


function  varargout  =  inputs(varargin) 

%  INPUTS  M-file  for  inputs. fig 

%  INPUTS,  by  itself,  creates  a  new  INPUTS  or  raises  the  existing 
%  singleton*. 

% 

%  H  =  INPUTS  returns  the  handle  to  a  new  INPUTS  or  the  handle  to 
%  the  existing  singleton*. 

% 
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%  INPUTS('CALLBACK',hObject,eventData, handles,...)  calls  the  local 

%  function  named  CALLBACK  in  INPUTS. M  with  the  given  input  arguments. 

% 

%  INPUTS('Property','Value',.--)  creates  a  new  INPUTS  or  raises  the 
%  existing  singleton*.  Starting  from  the  left,  property  value  pairs  are 
%  applied  to  the  GUI  before  inputs  OpeningFunction  gets  called.  An 
%  unrecognized  property  name  or  invalid  value  makes  property  application 
%  stop.  All  inputs  are  passed  to  inputs  OpeningFcn  via  varargin. 

% 

%  *See  GUI  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only  one 
%  instance  to  run  (singleton)". 

% 

%  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

%  Copyright  2002-2003  The  MathWorks,  Inc. 

%  Edit  the  above  text  to  modify  the  response  to  help  inputs 

%  Last  Modified  by  GUIDE  v2.5  10-Jan-2006  15:31:22 

%  Begin  initialization  code  -  DO  NOT  EDIT 
guiSingleton  =  1 ; 

gui_State  =  struct(’gui_Name',  mfilename, ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn',  @inputs_OpeningFcn,  ... 

'gui  OutputFcn',  @inputs_OutputFcn,  ... 

'gui_LayoutFcn',  [],... 

'guiCallback',  []); 
if  nargin  &&  ischar(varargin{l}) 

gui_State.gui_Callback  =  str2func(varargin{  1 }); 
end 

if  nargout 

[varargoutj  I  :nargout[]  =  gui_mainfcn(gui_State,  varargin):}); 
else 

gui_mainfcn(gui_State,  varargin  { : } ); 
end 

%  End  initialization  code  -  DO  NOT  EDIT 

%  —  Executes  just  before  inputs  is  made  visible. 

function  inputs_OpeningFcn(hObject,  eventdata,  handles,  varargin) 

%  This  function  has  no  output  args,  see  OutputFcn. 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  inputs  (see  VARARGIN) 
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%  Choose  default  command  line  output  for  inputs 
handles,  output  =  hObject; 

%  Update  handles  structure 
guidata(hObject,  handles); 

%  UIWAIT  makes  inputs  wait  for  user  response  (see  UIRESUME) 
%  uiwait(handles. figure  1); 


%  —  Outputs  from  this  function  are  returned  to  the  command  line, 
function  varargout  =  inputs_OutputFcn(hObject,  eventdata,  handles) 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Get  default  command  line  output  from  handles  structure 
varargout{l}  =  handles,  output; 

%  —  Executes  on  selection  change  in  popupmenul . 
function  popupmenul  _Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  popupmenul  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Convert  user  selection  of  coordinate  system  to  value  to  be  stored  as 
%  var.mat  (vehicle  pointing  system=  1 ,  ECI  system=2, 

%  Classic  orbital  elements  =3 
str=get(hObj  ect, 'String'); 
val=get(hObject, ’Value’); 
switch  sir  { val } ; 

case’Vehicle  Pointing’ 
var=l; 

case’Earth  Centered  Inertial’ 
var=2; 

case’Classic  Orbital  Elements' 
var=3; 
end 

save('var’); 

%  Hints:  contents  =  get(hObj ect, 'String')  returns  popupmenul  contents  as  cell  array 
%  contents  {get(hObj  ect, 'Value')}  returns  selected  item  from  popupmenul 
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%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  popupmenu  l_CreateFcn(hObj  ect,  eventdata,  handles) 

%  hObject  handle  to  popupmenul  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  popupmenu  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColorVwhite'); 
end 

%  —  Executes  on  selection  change  in  popupmenu2. 
function  popupmenu2_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  popupmenu2  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Convert  user  selection  of  atmosphere  shape  to  value  to  be  stored  as 
%  shape. mat  (spherical=l,  ellipsoid=2) 
str=get(hObject, 'String'); 
val=get(hObj  ect,’ V  alue’); 
switch  str  { val } ; 
case’Sphericaf 
shape=l; 

case’WGS-84  Ellipsoid’ 
shape=2; 
end 

save('shape’); 

%  Hints:  contents  =  get(hObj ect, ’String’)  returns  popupmenu2  contents  as  cell  array 
%  contents {get(hObject, ’Value’)}  returns  selected  item  from  popupmenu2 


%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  popupmenu2_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  popupmenu2  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  popupmenu  controls  usually  have  a  white  background  on  Windows. 
%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,  defaultUicontrolBackgroundColor’)) 
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set(hObject,’BackgroundColor’,'white'); 

end 

%  —  Executes  on  button  press  in  pushbutton  1. 

function  pushbutton  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  pushbutton  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Close  GUI  when  Run  is  selected 
close('inputs') 


Input  Function  for  Vehicle  Pointing  System 

The  file  vpoint.m  sets  up  a  graphical  interface  to  allow  the  user  to  input  the 
vehicle  parameters  for  mass,  reference  surface  area,  coefficient  of  lift,  and  coefficient  of 
drag.  It  also  allows  the  user  to  specify  the  initial  conditions  for  the  equations  of  motion: 
radius,  relative  flight-path  angle,  relative  velocity,  geocentric  latitude,  longitude,  relative 
heading,  roll  angle.  The  GUI  display  is  shown  in  Figure  (53)  followed  by  the  MATLAB 
code  for  v point. m. 
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Figure  61.  GUI  to  Input  Parameters  in  Vehicle  Pointing  System 


function  varargout  =  vpoint(varargin) 

%  VPOINT  M-file  for  vpoint.fig 

%  VPOINT,  by  itself,  creates  a  new  VPOINT  or  raises  the  existing 
%  singleton*. 

% 

%  H  =  VPOINT  returns  the  handle  to  a  new  VPOINT  or  the  handle  to 
%  the  existing  singleton*. 

% 

%  VPOINT('CALLBACK',hObject,eventData, handles,...)  calls  the  local 

%  function  named  CALLBACK  in  VPOINT.M  with  the  given  input  arguments. 

% 

%  VPOINT ('Property', 'Value',.  ••)  creates  a  new  VPOINT  or  raises  the 

%  existing  singleton*.  Starting  from  the  left,  property  value  pairs  are 
%  applied  to  the  GUI  before  vpoint  OpeningFunction  gets  called.  An 
%  unrecognized  property  name  or  invalid  value  makes  property  application 
%  stop.  All  inputs  are  passed  to  vpoint  OpeningFcn  via  varargin. 
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% 

%  *See  GUI  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only  one 
%  instance  to  run  (singleton)". 

% 

%  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

%  Copyright  2002-2003  The  MathWorks,  Inc. 

%  Edit  the  above  text  to  modify  the  response  to  help  vpoint 

%  Last  Modified  by  GUIDE  v2.5  10-Jan-2006  16:37:44 

%  Begin  initialization  code  -  DO  NOT  EDIT 
guiSingleton  =  1; 

gui_State  =  struct(’gui_Name',  mfilename, ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn',  @vpoint_OpeningFcn,  ... 

'gui  OutputFcn',  @vpoint_OutputFcn,  ... 

'gui_LayoutFcn',  [],... 

'gui_Callback',  []); 
if  nargin  &&  ischar(varargin{  1 }) 

gui_State.gui_Callback  =  str2func(varargin{  1 }); 
end 

if  nargout 

[varargout{l:nargout}]  =  gui_mainfcn(gui_State,  varargin{:}); 
else 

gui_mainfcn(gui_State,  varargin  { : } ); 
end 

%  End  initialization  code  -  DO  NOT  EDIT 

%  —  Executes  just  before  vpoint  is  made  visible. 

function  vpoint_OpeningFcn(hObject,  eventdata,  handles,  varargin) 

%  This  function  has  no  output  args,  see  OutputFcn. 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  vpoint  (see  VARARGIN) 

%  Choose  default  command  line  output  for  vpoint 
handles,  output  =  hObject; 

%  Update  handles  structure 
guidata(hObject,  handles); 

%  UIWAIT  makes  vpoint  wait  for  user  response  (see  UIRESUME) 
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%  uiwait(handles. figure  1); 

%  —  Outputs  from  this  function  are  returned  to  the  command  line, 
function  varargout  =  vpoint_OutputFcn(hObject,  eventdata,  handles) 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Get  default  command  line  output  from  handles  structure 
varargout)  1}  =  handles,  output; 

function  editl_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  as  a  double 

%  Convert  value  for  mass  to  matrix  and  save  to  m.rnat 

m=str2double(get(hObject, 'String')); 

save(’m') 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,'defaultUicontrolBackgroundColor’)) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  1  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  1  as  a  double 

%  Convert  value  for  coefficient  of  lift  to  matrix  and  save  to  Cl.mat 
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Cl=str2double(get(hObject, ’String')); 
save('Cl’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl  l_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0 , '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColorVwhite'); 
end 

function  editl2_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  editl2  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  12  as  a  double 

%  Convert  value  for  coefficient  of  drag  to  matrix  and  save  to  Cd.mat 

Cd=str2double(get(hObject, ’String’)); 

save(’Cd’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl2_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,’defaultUicontrolB  ackgroundColor’)) 
set(hObject,’BackgroundColor’, ’white’); 
end 

function  editl3_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 
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%  Hints:  get(hObject, 'String')  returns  contents  of  edit  13  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  13  as  a  double 

%  Convert  value  for  surface  area  to  matrix  and  save  to  S.mat 

S=str2double(get(hObject, 'String')); 

save('S'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl3_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0, defaultUicontrolBackgroundColor’)) 
set(hObject,'BackgroundColor','white'); 
end 

function  editl4_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  14  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATFAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  14  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  14  as  a  double 

%  Convert  value  for  radius  to  matrix  and  save  to  ri.mat 

ri=str2double(get(hObject, ’String’)); 

save('ri’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl4_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  14  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATFAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0,'defaultUieontrolBackgroundColor’)) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  editl5_Callback(hObject,  eventdata,  handles) 
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%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  15  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  15  as  a  double 

%  Convert  value  for  relative  flight-path  angle  to  matrix  and  save  to 
%  gammar.mat 

gammar=str2double(get(hObject, 'String')); 
save(’gammar'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl5_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,'defaultUieontrolBaekgroundColor')) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  editl6_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  16  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  16  as  a  double 

%  Convert  value  for  relative  velocity  to  matrix  and  save  to  vr.mat 

vr=str2double(get(hObject, 'String')); 

save(’vr'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl6_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER. 
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if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,  defaulfUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor','white'); 
end 

function  editl7_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  17  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  17  as  a  double 

%  Convert  value  for  geocentric  latitude  to  matrix  and  save  to  phii.mat 

phii=str2double(get(hObject, 'String’)); 

save(’phii’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl7_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,'defaultUicontrolBackgroundColor’)) 
set(hObject,’BackgroundColor’, 'white'); 
end 

function  editl8_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  18  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  1 8  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  18  as  a  double 

%  Convert  value  for  longitude  to  matrix  and  save  to  thetai.mat 

thetai=str2double(get(hObject, ’String’)); 

save(’thetai’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl8_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  18  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
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%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0, defaulfUicontrolBackgroundColor’)) 
set(hObject,’BackgroundColor','white’); 
end 

function  editl9_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl9  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  19  as  a  double 

%  Convert  value  for  relative  heading  to  matrix  and  save  to  psir.mat 

psir=str2double(get(hObject, 'String')); 

save('psir’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl9_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0, defaultUicontrolBackgroundColor’)) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  edit20_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit20  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit20  as  a  double 

%  Convert  value  for  bank  angle  to  matrix  and  save  to  sigma.mat 

sigma=str2double(get(hObject, ’String’)); 

save('sigma’); 
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%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit20_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor','white'); 
end 

%  —  Executes  on  button  press  in  pushbutton  1. 

function  pushbutton  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  pushbutton  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Close  GUI  when  RUN  pushbutton  is  pressed 
close('vpoinf); 


Input  Function  for  ECI  System 

The  file  eci.m  sets  up  a  graphical  interface  to  allow  the  user  to  input  the  vehicle 
parameters  for  mass,  reference  surface  area,  coefficient  of  lift,  and  coefficient  of  drag.  It 
also  allows  the  user  to  specify  the  initial  conditions  in  the  ECI  coordinate  system:  inertial 
position,  inertial  velocity,  roll  angle  and  Greenwich  sidereal  time.  The  GUI  display  is 
shown  in  Figure  (54)  followed  by  the  MATLAB  code  for  eci.m. 
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Figure  62.  GUI  to  Input  Parameters  in  ECI  System 


function  varargout  =  eci(varargin) 

%  ECI  M-file  for  eci.fig 

%  ECI,  by  itself,  creates  a  new  ECI  or  raises  the  existing 
%  singleton*. 

% 

%  H  =  ECI  returns  the  handle  to  a  new  ECI  or  the  handle  to 
%  the  existing  singleton*. 

% 

%  ECI('CALLBACK',hObject,eventData, handles,...)  calls  the  local 
%  function  named  CALLBACK  in  ECI.M  with  the  given  input  arguments. 
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% 

%  ECI('Property','Value',...)  creates  a  new  ECI  or  raises  the 
%  existing  singleton*.  Starting  from  the  left,  property  value  pairs  are 
%  applied  to  the  GUI  before  eciOpeningFunction  gets  called.  An 
%  unrecognized  property  name  or  invalid  value  makes  property  application 
%  stop.  All  inputs  are  passed  to  eci  OpeningFcn  via  varargin. 

% 

%  *See  GUI  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only  one 
%  instance  to  run  (singleton)". 

% 

%  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

%  Copyright  2002-2003  The  MathWorks,  Inc. 

%  Edit  the  above  text  to  modify  the  response  to  help  eci 

%  Last  Modified  by  GUIDE  v2.5  10-Jan-2006  18:1 1:33 

%  Begin  initialization  code  -  DO  NOT  EDIT 
guiSingleton  =  1 ; 

gui_State  =  struct('gui_Name',  mfilename, ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn',  @eci_OpeningFcn,  ... 

'gui  OutputFcn',  @eci_OutputFcn,  ... 

’gui_LayoutFcn',  [],... 

'guiCallback',  []); 
if  nargin  &&  ischar(varargin{  1 }) 

gui_State.gui_Callback  =  str2func(varargin{  1 }); 
end 

if  nargout 

[varargout{l:nargout}]  =  gui_mainfcn(gui_State,  varargin):}); 
else 

gui_mainfcn(gui_State,  varargin  { : } ); 
end 

%  End  initialization  code  -  DO  NOT  EDIT 

%  —  Executes  just  before  eci  is  made  visible. 

function  eci_OpeningFcn(hObject,  eventdata,  handles,  varargin) 

%  This  function  has  no  output  args,  see  OutputFcn. 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  eci  (see  VARARGIN) 

%  Choose  default  command  line  output  for  eci 
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handles,  output  =  hObject; 

%  Update  handles  structure 
guidata(hObject,  handles); 

%  UIWAIT  makes  eci  wait  for  user  response  (see  UIRESUME) 

%  uiwait(handles.  figure  1); 

%  —  Outputs  from  this  function  are  returned  to  the  command  line, 
function  varargout  =  eci_OutputFcn(hObject,  eventdata,  handles) 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Get  default  command  line  output  from  handles  structure 
varargout{l}  =  handles,  output; 

function  editl_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  as  a  double 

%  Convert  value  for  mass  to  matrix  and  save  as  m.mat 

m=str2double(get(hObject, 'String')); 

save('m') 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,  defaultUicontrolBackgroundColor’)) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
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%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  1  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  1  as  a  double 

%  Convert  value  for  coefficient  of  lift  to  matrix  and  save  as  Cl. mat 

Cl=str2double(get(hObject, 'String')); 

save('Cl’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl  l_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0,'defaultUicontrolBackgroundColor')) 
set(hObject,’BackgroundColor’,'white’); 
end 

function  editl2_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  12  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  12  as  a  double 

%  Convert  value  for  coefficient  of  drag  to  matrix  and  save  as  Cd.mat 

Cd=str2double(get(hObject, 'String')); 

save('Cd’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl2_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0,’defaultUicontrolBackgroundColor')) 
set(hObject,’BackgroundColor’,'white'); 
end 
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function  editl3_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  13  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  13  as  a  double 

%  Convert  value  for  reference  surface  area  to  matrix  and  save  as  S.mat 

S=str2double(get(hObject, 'String')); 

save('S’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl3_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl4_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  14  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  14  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  14  as  a  double 

%  Convert  value  for  ECI  x  position  to  matrix  and  save  as  x.mat 

x=str2double(get(hObject, 'String')); 

save('x'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl4_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  14  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER. 
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if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,  defaulfUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor','white'); 
end 

function  edit  15_Callback(hObj  ect,  eventdata,  handles) 

%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  15  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  15  as  a  double 

%  Convert  value  for  ECI  y  position  to  matrix  and  save  as  y.mat 

y=str2  double(get(hObj  ect,  ’  String’)) ; 

save('y'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit  15_CreateFcn(hObj  ect,  eventdata,  handles) 

%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,'defaultUicontrolBaekgroundColor’)) 
set(hObject,’BackgroundColor’, ’white’); 
end 

function  edit  16_Callback(hObj  ect,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  1 6  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  16  as  a  double 

%  Convert  value  for  ECI  z  position  to  matrix  and  save  as  z.mat 

z=str2double(get(hObject, 'String’)); 

save(’z'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit  16_CreateFcn(hObj  ect,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
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%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 


%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0, defaulfUicontrolBackgroundColor’)) 
set(hObject,’BackgroundColor','white’); 
end 

function  editl7_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  1 7  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  17  as  a  double 

%  Convert  value  for  ECI  x  velocity  to  matrix  and  save  as  vx.mat 

vx=str2double(get(hObject, 'String')); 

save('vx'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl7_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,defaultUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  editl8_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  18  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  1 8  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  18  as  a  double 

%  Convert  value  for  ECI  y  velocity  to  matrix  and  save  as  vy.mat 

vy=str2double(get(hObject, ’String’)); 

save('vy’); 
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%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl8_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  18  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColorVwhite'); 
end 

function  editl9_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  19  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  19  as  a  double 

%  Convert  value  for  ECI  z  velocity  to  matrix  and  save  as  vz.mat 

vz=str2double(get(hObject, ’String’)); 

save('vz'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl9_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,'defaultUieontrolBaekgroundColor’)) 
set(hObject,’BackgroundColor', ’white’); 
end 

function  edit20_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit20  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit20  as  a  double 
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%  Convert  value  for  bank  angle  to  matrix  and  save  as  sigma.mat 

sigma=str2double(get(hObject, 'String')); 

save('sigma'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit20_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,'BackgroundColor','white’); 
end 

function  edit21_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit21  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit21  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit21  as  a  double 

%  Convert  value  for  Greenwich  sidereal  hour  to  matrix  and  save  as  hour.mat 

hour=str2double(get(hObject, 'String')); 

save('hour'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit21_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit21  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,  defaultUicontrolBackgroundColor’)) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  edit22_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit22  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 
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%  Hints:  get(hObject, 'String')  returns  contents  of  edit22  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit22  as  a  double 

%  Convert  value  for  Greenwich  sidereal  minute  to  matrix  and  save  as  min.mat 

min=str2double(get(hObject, 'String')); 

save('min’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit22_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit22  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,'BackgroundColor','white'); 
end 

function  edit23_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit23  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit23  as  a  double 

%  Convert  value  for  Greenwich  sidereal  second  to  matrix  and  save  as  sec. mat 

sec=str2double(get(hObject, ’String’)); 

save(’sec’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit23_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 
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%  —  Executes  on  button  press  in  pushbutton  1. 

function  pushbutton  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  pushbutton  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Close  GUI  when  RUN  pushbutton  is  pressed 
close('eci'); 


Input  Function  for  Classic  Orbital  Elements 

The  file  orb.m  sets  up  a  graphical  interface  to  allow  the  user  to  input  the  vehicle 
parameters  for  mass,  reference  surface  area,  coefficient  of  lift,  and  coefficient  of  drag.  It 
also  allows  the  user  to  specify  the  initial  conditions  from  the  classic  orbital  elements: 
semi-major  axis,  true  anomaly,  inclination,  eccentricity,  argument  of  perigee,  right 
ascension  of  the  ascending  node,  roll  angle  and  Greenwich  sidereal  time.  The  GUI 
display  is  shown  in  Figure  (54)  followed  by  the  MATLAB  code  for  orb.m. 
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Figure  63.  GUI  to  Input  Parameters  from  Classic  Orbital  Elements 


function  varargout  =  orb(varargin) 

%  ORB  M-file  for  orb. fig 

%  ORB,  by  itself,  creates  a  new  ORB  or  raises  the  existing 
%  singleton*. 

% 

%  H  =  ORB  returns  the  handle  to  a  new  ORB  or  the  handle  to 
%  the  existing  singleton*. 

% 
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%  ORB('CALLBACK’,hObject,eventData, handles,...)  calls  the  local 

%  function  named  CALLBACK  in  ORB.M  with  the  given  input  arguments. 

% 

%  ORB('Property',' Value',. ••)  creates  a  new  ORB  or  raises  the 

%  existing  singleton*.  Starting  from  the  left,  property  value  pairs  are 
%  applied  to  the  GUI  before  orb  OpeningFunction  gets  called.  An 
%  unrecognized  property  name  or  invalid  value  makes  property  application 
%  stop.  All  inputs  are  passed  to  orb  OpeningFcn  via  varargin. 

% 

%  *See  GUI  Options  on  GUIDE'S  Tools  menu.  Choose  "GUI  allows  only  one 
%  instance  to  run  (singleton)". 

% 

%  See  also:  GUIDE,  GUIDATA,  GUIHANDLES 

%  Copyright  2002-2003  The  MathWorks,  Inc. 

%  Edit  the  above  text  to  modify  the  response  to  help  orb 

%  Last  Modified  by  GUIDE  v2.5  10-Jan-2006  18:07:37 

%  Begin  initialization  code  -  DO  NOT  EDIT 
guiSingleton  =  1 ; 

gui_State  =  struct(’gui_Name',  mfilename, ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn',  @orb_OpeningFcn,  ... 

'gui  OutputFcn',  @orb_OutputFcn,  ... 

'gui_LayoutFcn',  [],... 

'guiCallback',  []); 
if  nargin  &&  ischar(varargin{l}) 

gui_State.gui_Callback  =  str2func(varargin{  1 }); 
end 

if  nargout 

[varargoutj  I  margout}]  =  gui_mainfcn(gui_State,  varargin):}); 
else 

gui_mainfcn(gui_State,  varargin  { : } ); 
end 

%  End  initialization  code  -  DO  NOT  EDIT 

%  —  Executes  just  before  orb  is  made  visible. 

function  orb_OpeningFcn(hObject,  eventdata,  handles,  varargin) 

%  This  function  has  no  output  args,  see  OutputFcn. 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATL AB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  orb  (see  VARARGIN) 
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%  Choose  default  command  line  output  for  orb 
handles,  output  =  hObject; 

%  Update  handles  structure 
guidata(hObject,  handles); 

%  UIWAIT  makes  orb  wait  for  user  response  (see  UIRESUME) 

%  uiwait(handles. figure  1); 

%  —  Outputs  from  this  function  are  returned  to  the  command  line, 
function  varargout  =  orb_OutputFcn(hObject,  eventdata,  handles) 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Get  default  command  line  output  from  handles  structure 
varargout)  1}  =  handles,  output; 

function  editl_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  as  a  double 

%  Convert  value  for  mass  to  matrix  and  save  as  m.mat 

m=str2double(get(hObject, 'String')); 

save('m') 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl  l_Callback(hObject,  eventdata,  handles) 
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%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  editl  1  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  editl  1  as  a  double 

%  Convert  value  for  coefficient  of  lift  to  matrix  and  save  as  Cl. mat 

Cl=str2double(get(hObject, 'String')); 

save('Cl'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl  l_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl2_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  12  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  12  as  a  double 

%  Convert  value  for  coefficient  of  drag  to  matrix  and  save  as  Cd.mat 

Cd=str2double(get(hObject, ’String’)); 

save(’Cd’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl2_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  12  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER. 

if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 

get(0,  defaulfUicontrolBackgroundColor')) 
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set(hObject,’BackgroundColor’,'white'); 

end 

function  editl3_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String')  returns  contents  of  edit  13  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  13  as  a  double 

%  Convert  value  for  reference  surface  area  to  matrix  and  save  as  S.mat 

S=str2double(get(hObject, 'String')); 

save(’S’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl3_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  13  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0, defaultUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor','white’); 
end 

function  editl5_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  1 5  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  15  as  a  double 

%  Convert  value  for  semi-major  axis  to  matrix  and  save  as  a.mat 

a=str2double(get(hObject, 'String')); 

save('a’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl5_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  15  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 
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%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,'BackgroundColorVwhite’); 
end 

function  editl6_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String')  returns  contents  of  edit  16  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  16  as  a  double 

%  Convert  value  for  true  anomaly  to  matrix  and  save  as  nu.mat 

nu=str2double(get(hObject, 'String')); 

save('nu’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl6_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  16  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,'defaultUieontrolBackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 

function  editl7_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit  1 7  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit  17  as  a  double 

%  Convert  value  for  inclination  to  matrix  and  save  as  inc.mat 

inc=str2double(get(hObject, 'String')); 

save('inc'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl7_CreateFcn(hObject,  eventdata,  handles) 
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%  hObject  handle  to  edit  17  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,'BackgroundColorVwhite'); 
end 

function  editl8_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  18  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  18  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  18  as  a  double 

%  Convert  value  for  argument  of  perigee  to  matrix  and  save  as  omega.mat 

omega=str2double(get(hObject, ’String’)); 

save(’omega’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl8_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  editl  8  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor’), 
get(0,'defaultUicontrolB  ackgroundColor’)) 
set(hObject,’BackgroundColor', ’white’); 
end 

function  editl9_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit  19  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit  19  as  a  double 

%  Convert  value  for  right  ascension  of  ascending  node  to  matrix  and  save  as 
%  node. mat 
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node=str2double(get(hObject, ’String’)); 
save('node’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  editl9_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit  19  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColorVwhite’); 
end 

function  edit20_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit20  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit20  as  a  double 

%  Convert  value  for  bank  angle  to  matrix  and  save  as  sigma.mat 

sigma=str2double(get(hObject, 'String')); 

save('sigma’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit20_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit20  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,’BackgroundColor'), 
get(0,’defaultUicontrolB  ackgroundColor’)) 
set(hObject,’BackgroundColor’, 'white'); 
end 

function  edit21_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit21  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 
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%  Hints:  get(hObject, 'String')  returns  contents  of  edit21  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit21  as  a  double 

%  Convert  value  for  Greenwich  sidereal  hour  to  matrix  and  save  as  hour.mat 

hour=str2double(get(hObject, 'String')); 

save('hour'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit21_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit21  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0, defaultUicontrolBackgroundColor’)) 
set(hObject,'BackgroundColor','white'); 
end 

function  edit22_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit22  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATEAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String’)  returns  contents  of  edit22  as  text 
%  str2double(get(hObject, ’String’))  returns  contents  of  edit22  as  a  double 

%  Convert  value  for  Greenwich  sidereal  minute  to  matrix  and  save  as  min. mat 

min=str2double(get(hObject, ’String’)); 

save('min'); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit22_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit22  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,'defaulfUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor', 'white'); 
end 

function  edit23_Callback(hObject,  eventdata,  handles) 


124 


%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, 'String')  returns  contents  of  edit23  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit23  as  a  double 

%  Convert  value  for  Greenwich  sidereal  second  to  matrix  and  save  as  sec. mat 

sec=str2double(get(hObject, 'String')); 

save('sec’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit23_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER, 
if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0 , ' '  defaultU  icontrolB  ackgroundColor')) 
set(hObject,’BackgroundColor', 'white'); 
end 


function  edit24_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hints:  get(hObject, ’String')  returns  contents  of  edit23  as  text 
%  str2double(get(hObject, 'String'))  returns  contents  of  edit23  as  a  double 

%  Convert  value  for  eccentricity  to  matrix  and  save  as  e.mat 

e=str2double(get(hObject, 'String')); 

save('e’); 

%  —  Executes  during  object  creation,  after  setting  all  properties, 
function  edit24_CreateFcn(hObject,  eventdata,  handles) 

%  hObject  handle  to  edit23  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

%  Hint:  edit  controls  usually  have  a  white  background  on  Windows. 

%  See  ISPC  and  COMPUTER. 
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if  ispc  &&  isequal(get(hObject,'BackgroundColor'), 
get(0,  defaultUicontrolBackgroundColor')) 
set(hObject,'BackgroundColor’,'white'); 
end 

%  —  Executes  on  button  press  in  pushbutton  1. 

function  pushbutton  l_Callback(hObject,  eventdata,  handles) 

%  hObject  handle  to  pushbutton  1  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Close  GUI  when  RUN  pushbutton  is  pressed 
close(’orb'); 


Ordinary  Differential  Equation  Integrators 

The  core  program  reentry. m  calls  one  of  two  functions  to  perform  the  integration 
of  the  equations  of  motion.  The  function  spherical,  m  is  called  when  the  atmosphere 
selected  by  the  user  is  spherical  and  the  function  ellipsoid. m  is  called  when  the  user 
selects  an  ellipsoidal  atmosphere.  In  both  cases  Equations  (36)  and  (37)  for  calculating 
drag  and  lift  are  embedded  in  the  equations  of  motion.  In  the  function  ellipsoid. m 
Equations  (44)  and  (45),  used  to  calculate  the  Earth’s  radius  at  a  given  geocentric 
latitude,  are  also  embedded  into  the  equations  of  motion. 


Code  for  Spherical  Atmosphere  Integrator 

%  This  integration  function  is  called  by  reentry. m  when  a  spherical  atmosphere  is 
%  selected 

function  dy=reentry(t,y,m,S, Cl, Cd, sigma) 

%  Earth’s  acceleration  from  gravity  at  equatorial  radius  (m/sA2) 
g0=9.798; 

%  Surface  density  of  atmosphere  at  sea  level  (kg/mA3) 
rhos=l. 54793; 

%  Radius  of  Earth  (m) 
re=6378137; 
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%  Inverse  of  the  scale  height  of  atmosphere  (1/m) 
beta=l/6935; 

%  WGS-84  rotation  rate  of  Earth  (rad/s) 
omegae=. 000072921 151467; 

%  Establishing  dy  vector  for  equations  of  motion 
dy=zeros(6,l); 

%  The  equations  of  motion 

%  dr/dt  (solved  for  radius) 
dy(l)=y(2)*sin(y(3)); 

%  dV/dt  (solved  for  relative  velocity) 

dy(2)=-((rhos*  exp(-beta*  (y(  1  )-re))  *  S  *Cd*y(2)A2)/2)/m-gO  *  (re/y(  1  ))A2  * . . 
sin(y(3))+y(l)*omegaeA2*cos(y(3))*(cos(y(5))*sin(y(3))-sin(y(5))... 
*sin(y(6))*cos(y(3))); 

%  dgamma/dt  (solved  for  relative  flight-path  angle) 
dy(3)=((rhos*  exp(-beta*  (y(  1  )-re))  *  S  *  Cl*y(2) A2)/2)  *  cos(sigma)/ (y(2)  *m) 
g0*(re/y(l))A2*cos(y(3))/y(2)+y(2)*cos(y(3))/y(l)+y(l)*omegaeA2*... 
cos(y(5))*(cos(y(5))*cos(y(3))+sin(y(5))*sin(y(6))*sin(y(3)))/y(2)... 

+2*omegae*cos(y(5))*cos(y(6)); 

%  dtheta/dt  (soved  for  longitude) 
dy(4)=y(2)*cos(y(3))*cos(y(6))/(y(l)*cos(y(5))); 

%  dphi/dt  (solved  for  latitude) 
dy(5)=y(2)*cos(y(3))*sin(y(6))/y(l); 

%  dpsi/dt  (solved  for  relative  heading) 

dy(6)=((rhos  *  exp(-beta*  (y(  1  )-re))  *  S  *  Cl*y(2)A2)/2)*  sin(sigma)/ (m* . . . 
cos(y(3))*y(2))-y(2)*cos(y(3))*cos(y(6))*tan(y(5))/y(l)+2s|somegae5|!... 
(sin(y(6))*cos(y(5))*tan(y(3))-sin(y(5)))-y(l)*omegaeA2*sin(y(5))*... 
cos(y(5))*cos(y(6))/(cos(y(3))5|!y(2)); 


Code  for  Ellipsoidal  Atmosphere  Integrator 

%  This  integration  function  is  called  by  reentry. m  when  an  ellipsoidal 

%  atmosphere  is  selected 

function  dy=reentry(t,y,m,S, Cl, Cd, sigma) 

%  Earth's  acceleration  due  to  gravity  at  equatorial  radius  (m/sA2) 
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g0=9.798; 


%  Earth's  atmosphere’s  surface  density  (kg/mA3) 
rhos=l. 54793; 

%  WGS-84  value  for  the  radius  of  the  Earth  at  the  equator  (m) 
re=6378137; 

%  WGS-84  value  for  the  radius  of  the  Earth  at  the  poles  (m) 
rp=6356752.3142; 

%  Ellipticity  of  the  Earth  (e) 
e=l-rp/re; 


%  Inverse  of  the  scale  height  for  Earth's  atmosphere  (1/m) 
beta=l/6935; 

%  WGS-84  value  for  Earth's  rotation  rate  (rad/s) 
omegae=. 000072921 151467; 

%  Establishing  dy  vector  for  equations  of  motion 
dy=zeros(6,l); 

%  The  equations  of  motion 

%  dr/dt  (solved  for  radius) 

dy(l)=y(2)*sin(y(3)); 

%  dV/dt  (solved  for  relative  velocity) 

dy(2)=-((rhos  *  exp(-beta*(y(  1  )-rp/(cos(y(5))  *((tan(y(5)))A2+. . . 

(l-e)A2)A.5)))*S*Cd*y(2)A2)/2)/m-g0*(re/y(l))A2*sin(y(3))+y(l)*... 

omegaeA2*cos(y(3))*(cos(y(5))*sin(y(3))-sin(y(5))*sin(y(6))*cos(y(3))); 

%  dgamma/dt  (solved  for  relative  flight-path  angle) 
dy(3)=((rhos*exp(-beta*(y(l)-rp/(cos(y(5))*((tan(y(5)))A2+... 

(l-e)A2)A.5)))*S*Cl*y(2)A2)/2)*cos(sigma)/(y(2)*m)-g0*(re/y(l))A2*... 
cos(y(3))/y(2)+y(2)5|!cos(y(3))/y(T)+y(l)5|!omegaeA2s|5cos(y(5))!,=  ... 
(cos(y(5))*cos(y(3))+sin(y(5))*sin(y(6))*sin(y(3)))/y(2)+2*omegae*... 
cos(y(5))*cos(y(6)); 

%  dtheta/dt  (solved  for  longitude) 
dy(4)=y(2)*cos(y(3))*cos(y(6))/(y(l)*cos(y(5))); 

%  dphi/dt  (solved  for  geocentric  latitude) 
dy(5)=y(2)*cos(y(3))5|!sin(y(6))/y(l); 
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%  dpsi/dt  (solved  for  relative  heading) 
dy(6)=((rhos*exp(-beta*(y(l)-rp/(cos(y(5))*((tan(y(5)))A2+... 

( 1  -e)A2)A .  5)))  *  S  *  Cl*y(2)A2)/2)  *  sin(sigma)/ (m*  cos(y(3))  *y(2))-y(2)* . . . 
cos(y(3))*cos(y(6))*tan(y(5))/y(l)+2*omegae*(sin(y(6))*cos(y(5))*... 
tan(y(3))-sin(y(5)))-y(l)*omegaeA2*sin(y(5))*cos(y(5))*cos(y(6))/... 
(cos(y(3))*y(2)); 


Functions  to  Terminate  Integration 

The  final  two  functions  in  the  MATLAB  software  package  are  used  to  halt  the 
integration  functions  when  the  entry  vehicle  reaches  the  surface  of  the  Earth.  The 
function  stops.m  is  called  by  the  core  program  reentry. m  when  a  spherical  atmosphere  is 
selected  by  the  user.  The  function  is  established  by  the  odeset  command  as  an  “Event” 
and  stops  the  integration  function  spherical,  m  when  the  radius  solved  for  in  the  first 
equation  of  motion  is  less  than  the  constant  radius  of  the  Earth.  Likewise,  the  “Event” 
stope.m  is  implemented  when  an  ellipsoidal  atmosphere  is  selected,  stope.m  halts 
integration  when  the  radius  solved  for  in  the  first  equation  of  motion  is  less  than  the 
radius  of  the  Earth  found  in  Equations  (44)  and  (45)  as  a  function  of  geocentric  latitude. 


Termination  for  Spherical  Earth 

%stops.m 

%Events  function  to  halt  integration  at  radius  of  spherical  Earth 

function  [value, isterminal, direction]  =  stops(t,y,m,S, Cl, Cd, sigma) 

%  Detect  point  that  radius  is  less  than  radius  of  Earth 
if  y(l)  <  6378137; 

value  =  0; 
else 

value  =  1 ; 
end 

%  Halt  integration  when  value=0  is  detected 
isterminal  =  1 ; 
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%  Halt  for  all  value=0 
direction  =  0; 


Termination  for  Ellipsoidal  Earth 

%stope.m 

%Events  function  to  halt  integration  at  radius  of  ellipsoidal  Earth 

function  [value, isterminal, direction]  =  stope(t, y,m,S, Cl, Cd, sigma) 

%  Detect  point  that  radius  is  less  than  radius  of  Earth 
ify(l)  <  6356752.3 142/(cos(y(5))*((tan(y(5)))A2+... 

(6356752.3 142/6378 137)A2)A.5); 
value  =  0; 
else 

value  =  1 ; 
end 

%  Halt  integration  when  value=0  is  detected 
isterminal  =  1 ; 

%  Halt  for  all  value=0 
direction  =  0; 
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